The ANI-MBIS-q model was trained on precalculated MBIS charges from structures in the ANI-2x dataset, which correspond to electronic densities converged at the ωB97X/def2-TZVPP level of theory. The basis set used is slightly larger than the one used to train the original ANI-2x model (6-31G*), we note however that MBIS charges have been shown to be largely independent of basis set size. In this work, we present and showcase the implementation of a ML/MM scheme that takes advantage of the capability of the ANI-MBIS-q model to predict accurate atomic charges and charge-gradients, together with that of the ANI-2x model of producing fast and accurate energies and forces at the ωB97X/6-31G* level of theory, by combining both together into a single model, ANI-MBIS, and embedding it in the context of a previously thoroughly validated ML/MM coupling approach. In addition, we introduce an interface that allows us to use ANI-style neural networks implemented in the TorchANI 15 framework, as calculators for the QM region in QM/MM simulations performed in SANDER, the molecular dynamics (MD) engine provided in AmberTools, one of the leading software-suites for biomolecular simulations. This implementation, which is planned for inclusion in the next AmberTools release, is efficient, in that it avoids the cost that would incur from the naive approach of calling into Python code for every model evaluation. It also requires no familiarity with the internals of PyTorch, or any other machine learning framework from prospective users, and can be seamless called through the SANDER engine. We finally note that, the ANI models being general-purpose models that perform well for a wide range of systems, our scheme can be used without the need to train models from scratch. This article is organized as follows: We first discuss the theory and implementation details of our ML/MM scheme. Next, we evaluate the performance and accuracy of our implementation for obtaining accurate MD simulations by the calculation of solvation profiles, vibrational spectra through the vibrational density of states (VDOS), and torsional free-energy profiles of small molecules in aqueous solution. We then assess our ML/MM approach's capability to generate a physically accurate description of protein-ligand interactions. Finally, we analyze the performance of our implementation and offer perspectives for future improvements.