Efficient Calculation of Electrostatic Energies for Large-Scale Nonadiabatic Molecular Dynamics in a Site Basis

Nonadiabatic molecular dynamics simulation of charge and exciton transport in molecular materials and biological systems are often carried out in a (quasi-)diabatic or site basis. Such simulations require the calculation of the electrostatic site energy of all possible charge or excited states of the system at each molecular dynamics step, which quickly becomes computationally prohibitive when Ewald summation is used. By combining the damped shifted force real space electrostatic summation method with a suitable addition-subtraction scheme, we show that the calculation of electrostatic energy and forces for N_mol site energies can be carried out at a small and system size independent overhead compared to the calculation for a single site energy. This advance enables us to include full electrostatic interactions in nonadiabatic molecular dynamics simulations for charge and exciton transport. Applying our computational scheme to hole transport in crystalline anthracene, we find that upon inclusion of electrostatic site energy fluctuations (also sometimes termed diagonal electrostatic disorder) the inverse participation ratio measuring hole delocalization decreases from ∼5 to ∼4 concomitant with a decrease in the hole mobility by about 9% along the b-crystallographic direction and by 30% along the a-direction. Accounting for electrostatics improves the agreement with experimental time-of-flight mobilities and mobility anisotropy, but it does not alter the charge transport mechanism, transient delocalization. Our work confirms that omission of electrostatic site energy disorder is a reasonable approximation for acenes, yet electrostatics is required to obtain near-quantitative agreement with experiment, even for apolar systems.