This work proposes to simulate the direct photoionization process of small aromatic clusters of benzene, naphthalene, and anthracene, using the combination of a tailored model potential and a valence bond Hamiltonian. Classical molecular dynamics simulations are carried out on the potential energy surface of the neutral cluster, at energies meant to transcribe quantum mechanical effects such as intermolecular vibrational zero-point energy and the classically forbidden area of the vibrational ground-state wave function. The statistics of the ionization potential values obtained by vertical ionization of neutral configurations originating from these simulations yield valuable information on the respective topologies of the neutral and ionic potential energy surfaces. Comparison of calculated values with experimental one-photon and two-photon ionization potential determinations and ion current curves is also performed.