The objective of this research is to properly calculate the association energy with molecular models of asphaltenes using Molecular Dynamics (MD) simulations. Asphaltene precipitation clogs pore throats of reservoir rocks, and significantly reduces the production rate. To better understand the asphaltene precipitation, a process-modeling methodology has been developed. The cubic-plus association equation of state (CPA) (Li & Firoozabadi, 2010) is effective in the modeling, and contains only one adjustable parameter, that is, association energy of asphaltene, however, we were not able to correctly evaluate the association energy by the currently available methods. Quantitative Molecular Representations (QMR) method provides systematic structural representations of asphaltenes based on analytical data of a crude oil, from which molecular models of asphaltene can be obtained. These models will be employed in this study for calculation of the association energy. Two different asphaltene models were generated for two different crude oils (Oil 'A' and Oil 'M'). The potential of mean force (PMF) is calculated with two models of asphaltene molecules in liquid benzene for two kinds of asphaltenes using umbrella sampling method. The association energy is calculated by the difference between values at the lowest bottom of a trough and at infinity of the PMF curve. As the calculation conditions, temperature and pressure are 300 K and 100 bar, respectively. The PMF profiles for asphaltenes from Oil 'A' and Oil 'M' were obtained, and the values of association energy were estimated as -8.39 kJ/mol and -8.91 kJ/mol, respectively. Two aromatic planes of asphaltene molecules are stacked with each other at the minimum energy point. This π-π stacking is one of the main interactions to form asphaltene nano-aggregation. The results of asphaltenes from Oil 'A' and Oil 'M' have a similar tendency; for example, the distance between the centers of mass of two asphaltenes at the minimum energy points factually equal. The data obtained in this work should be used in the thermodynamic model of phase behavior of asphaltene dealing with asphaltene association and precipitation. It will improve the prediction ability of reservoir simulators with cubic-plus association equation of state.