We report a study on the non-adiabatic molecular dynamics (NA-MD) of the charge transfer (CT) process in the boron subphtalocyanine chloride (SubPc)/fullerene (C60) interface using our newly implemented Libra-X software package, which is based on an interface of the Libra NA-MD library and the GAMESS electronic structure software. In particular, we address the following aspects of the simulation protocol: (a) the choice of the potential used to treat interatomic interactions and its effect on the structures of the complex and CT rates; (b) the choice of the electronic structure methodology used; and (c) the choice of the trajectory surface hopping (TSH) methodology used. From our analysis of the electronic structure, we suggest that the distortion of the SubPc conical structure affects orbital localization and that the "breathing" motion of SubPc drives the CT process in SubPc/C60. This study illustrates that the choice of the TSH methodology and electronic decoherence are crucial for the CT simulation. We extend our analysis of CT in SubPc/(C60)n models by increasing the number of C60 molecules up to n = 4. We find that the details of the interfacial SubPc/(C60)n geometry determine the CT rate. Finally, we find the computed CT timescale to be in the range of 2.2-5.0 ps, which is in agreement with the experimentally determined timescale in the order of magnitude of ∼10 ps. The developed open-source Libra-X package is freely available on the Internet at https://github.com/Quantum-Dynamics-Hub/Libra-X.