High-Throughput Docking Using Quantum Mechanical Scoring

Today high-throughput docking is one of the most commonly used computational tools in drug lead discovery. While there has been an impressive methodological improvement in docking accuracy, docking scoring still remains an open challenge. Most docking programs are rooted in classical molecular mechanics. However, to better characterize protein-ligand interactions, the use of a more accurate quantum mechanical (QM) description would be necessary. In this work, we introduce a QM-based docking scoring function for high-throughput docking and evaluate it on 10 protein systems belonging to diverse protein families, and with different binding site characteristics. Outstanding results were obtained, with our QM scoring function displaying much higher enrichment (screening power) than a traditional docking method. It is acknowledged that developments in quantum mechanics theory, algorithms and computer hardware throughout the upcoming years will allow semi-empirical (or low-cost) quantum mechanical methods to slowly replace force-field calculations. It is thus urgently needed to develop and validate novel quantum mechanical-based scoring functions for high-throughput docking toward more accurate methods for the identification and optimization of modulators of pharmaceutically relevant targets.


INTRODUCTION
The cost to bring a new drug to the market could be as high as 2.6 billion US dollars, and can take up to 15 years (DiMasi et al., 2016). For many years, both the identification and optimization of novel drug lead compounds were accomplished within the drug discovery process by the experimental high-throughput screening of large chemical libraries. In spite of multiple efforts to improve its performance, drug discovery remains a costly and time consuming technique (Phatak et al., 2009). However, for the last 25 years, theoretical developments, better computational algorithms, faster computing resources, and improved visualization tools enabled the routine use of computational methods to model and visualize protein-ligand (PL) interactions, calculate binding free energy to different degrees of accuracy, and in silico screen chemical libraries using ligand-based and structure-based approaches. Today, computational chemistry is firmly established as a valuable tool in any drug lead discovery endeavor, aimed at saving time, effort, resources, and reducing costs (Cavasotto and Orry, 2007;Jorgensen, 2009Jorgensen, , 2012Spyrakis and Cavasotto, 2015;Pagadala et al., 2017).
During the last three decades, molecular docking has been one of the most commonly used computational methods in drug lead discovery (for review, cf., Kitchen et al., 2004;Rognan, 2011;Ciancetta and Moro, 2015;Sotriffer, 2015;Spyrakis and Cavasotto, 2015;Sulimov et al., 2019b). The aim of proteinsmall-molecule docking is the characterization of the optimal binding modes (poses) of a molecule within the binding site, and an estimation of its binding free energy. In high-throughput docking (HTD), where the protein is usually considered rigid or with very few degrees of freedom, and thousands to millions of molecules from a chemical library are screened, the goal is to generate a sub-library enriched with potential ligands, which will be prioritized for further experimental evaluation. In HTD, two different stages can be distinguished: the assessment of the best binding mode(s) of each molecule of the library ("docking stage"), and, on each in silico generated protein-small-molecule complex, the calculation of a score reflective of the likelihood that the molecule will actually bind to the target ("scoring stage") (Cavasotto and Orry, 2007;Guedes et al., 2018). In the docking stage, the docking energy (DE) is used to select, for each molecule, the lowest-energy pose(s) from a large amount of conformations generated, while the docking score (DS) is generally calculated as a fast approximation to the binding free energy ( G binding ), and depends on several factors, such as the energy representation of the system, the model used to represent the aqueous environment and the consideration of explicit water molecules within the active site (Cozzini et al., 2006;Amadasi et al., 2008), and the degree of consideration of receptor flexibility (Cavasotto and Singh, 2008;Spyrakis et al., 2011;Spyrakis and Cavasotto, 2015). Thus, DE discriminates among poses of the same molecule, while the DS characterizes each molecule of the docked chemical library and is used to rank them according to the likelihood of binding. Many docking programs, however, use a single function as DE and DS.
It should be stressed that one of the main advantages of docking is that in silico generated poses usually serve as the starting point for in silico ligand optimization, using for example molecular dynamics-based calculation of binding free energies, such as Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) and MM-Generalized Born Surface Area (MM-GBSA) methods (Kerrigan, 2013;Reddy et al., 2014;Genheden and Ryde, 2015;Sun et al., 2018;Wang et al., 2019).
While docking accuracy depends on the program, it is acknowledged that most of them are usually successful in identifying the correct pose (RMSD < 2 Å) with respect to the native structure (Warren et al., 2006;Wang et al., 2016). Moreover, an extensive recent benchmark of the Comparative Assessment of Scoring Functions (CASF) (Su et al., 2019) highlighted that docking programs display a better performance in terms of docking accuracy than in any of these three scoringrelated metrics: correlation with experimental binding data (scoring power), ranking of ligands by their binding affinity data provided their correct poses are known (ranking power), and identification of actual ligands from a sub-library of top-ranking small-molecules (screening power). This was in agreement with other works (Cavasotto and Abagyan, 2004;Slater and Kontoyianni, 2019).
Most docking developments have been mainly rooted in molecular mechanics (MM) force-fields (FF). However, to better characterize protein-ligand interactions, at least in some cases, the use of a quantum mechanical (QM) description would be necessary . The QM formulation is theoretically exact, as in principle, it accounts for all contributions to the energy (including terms or effects usually missing in FFs, such as electronic polarization, charge transfer, halogen bonding, and covalent-bond formation). Moreover, the QM framework is general across the chemical space so that all elements and interactions can be considered on equal footing, thus avoiding MM parameterizations.
Following the pioneering work of Merz (2004, 2005) where a QM-based score was used to discriminate ligand from decoy poses, there have been recently some applications of QM methods in docking, mainly aiming for accurate ligand binding mode assessment (for a survey of recent related works cf., Mucs and Bryce, 2013;Cavasotto et al., 2018;Aucar and Cavasotto, 2020). In a significant step forward, Pecina et al. obtained impressive results on the discrimination of native from decoy docking poses on four challenging systems (Pecina et al., 2016) using a docking energy function (Lepšík et al., 2013) based on the semi-empirical QM PM6 Hamiltonian (Stewart, 2007) supplemented with the D3H4X correction for dispersion, hydrogen-and halogenbonding interactions (Rezáč and Hobza, 2012). In a follow-up contribution (Pecina et al., 2017), an even superior performance was achieved for accurate pose assessment using a self-consistentcharge density-functional tight-binding method (SCC-DFTB) formulation coupled with D3H4 corrections for dispersion and hydrogen-bond interactions, though at a higher computational cost. This docking energy score function was further used to obtain a reliable ranking on 10 inhibitors binding to carbonic anhydrase II (CAII) (Pecina et al., 2018).
However, the development of QM-based docking scoring functions aiming at the ranking of molecules within HTD (screening power) has progressed at a significantly slower pace. Only very recently, a QM-based approach was presented displaying a very good performance on discriminating ligands and decoys on a single system (heat shock protein 90, HSP90) (Eyrilmez et al., 2019). In fact, the development of fast yet accurate docking scoring functions still constitutes an area of active research (Cavasotto, 2012;Guedes et al., 2018). Moreover, the blind challenges ran by the Drug Design Data Resource (D3R) for ligand-pose and affinity prediction in 2015 (Gathiaka et al., 2016), 2016 (Gaieb et al., 2018), and 2018 (Gaieb et al., 2019), have shown the importance of method development and benchmarking in pose prediction and binding affinity ranking of ligands.
In this work, we introduce a QM-based docking scoring function and evaluate it in terms of ligand enrichment on 10 protein systems belonging to diverse protein families in terms of different binding site characteristics, the presence of co-factors and water molecules, and the enrichment factors computed with a standard HTD method. Excellent results were obtained by displaying our QM-based scoring function a much higher enrichment (screening power) than a traditional docking method. We stress that our goal is to present and to validate an initial straightforward approach, which could serve as a starting point for further developments and improvement. A wider and extensive benchmarking on more systems and a systematic comparison with most of the standard docking programs, and the assessment of the optimal combination of the different components of our approach (QM formalism and continuum solvent model, energy minimization strategies, use of single or multiple docking poses for scoring, and entropy contribution) are considerations of their importance. However, they exceed the purpose of our work and will be published in due course.
Assuming a continuous development in QM theory, algorithms and computer hardware, it is likely that semiempirical methods [or low-cost Density Functional Theory (DFT) methods] will replace FF over the next 25 years (Grimme and Schreiner, 2018). Therefore, it is absolutely justified and there is an urgent need to start developing the next generation of QMbased scoring functions for HTD toward better methods for the identification of small-molecule modulators of pharmaceutically relevant targets.

Protein Systems Preparation
The following targets were downloaded from the PDB (cf. Table 1) Each target was prepared using ICM software (MolSoft, San Diego, CA, 2019;  in a similar fashion as in earlier works (Phatak et al., 2010). Succinctly, hydrogen atoms were added, followed by a local energy minimization of the complete system, and polar and water hydrogen positions were determined by optimizing the hydrogen bonding network within the torsional coordinates space. All Asp and Glu residues were assigned a −1 charge, and all Arg and Lys residues were assigned a +1 charge. Histidine tautomers were chosen according to their corresponding hydrogen bonding pattern. For docking with AutoDock Vina (Trott and Olson, 2010), the systems were pre-processed with AutoDock Tools (Morris et al., 2009).

Docking Library Preparation
For each target, the docking libraries were built by merging a set of ligands and a set of decoys, where the latter had similar physico-chemical properties to the ligands, but dissimilar 2-D topology. This has been shown to be necessary to ensure unbiased results when benchmarking docking programs (Huang et al., 2006;Gatica and Cavasotto, 2012). Ligands and decoys were extracted from the Directory of Useful Decoys (DUD, Huang et al., 2006), the NRLiSt binding data base for nuclear receptors (Lagarde et al., 2014),

Protein-Molecule Complex Generation, Structural Relaxation, and Unbound Protein and Ligand States Characterization
Protein-molecule complexes for QM-scoring were generated using the ICM docking module, keeping for each molecule its lowest DE conformation (docking RMSD values of native ligands are shown in Table 2). These protein-molecule complexes were also relaxed through cycles of local energy minimization in ICM according to the following procedure: (i) For each protein, the collected dihedral angles of amino-acids within 4 Å of any docked ligand of the corresponding chemical library were considered free; (ii) For each protein-molecule complex, five cycles of local energy minimization were performed restraining the heavy atoms with a harmonic potential with respect to their initial conformation; in each cycle the weight of this added potential was reduced in the following way: 50, 10, 5, 1, and 0 kcal/mol (no restraint). During this local energy minimization, the protein system was optimized in the torsional space , and the small-molecule in the Cartesian space.
To generate the unbound states, local energy minimization was performed on both protein and small-molecule in isolation from the crystallographic structure and the docked conformation, respectively.

System Cutout
For each target, a reduced-system was defined by first listing all the amino-acids within 8 Å of any docked molecule with ICM (only heavy atoms were considered in this threshold). Then, upon visual inspection, other amino-acids were eventually added to the list in order to avoid intra-helix or intra-β-sheet fragmentation, or loop fragments with just one amino-acid. A reduced-system was then built by deleting from the structure all amino-acids not included in the list, capping the N-and C-terminal of each fragment with hydrogens.

Entropy Calculation
Binding small-molecule conformational entropy was estimated as where it is assumed that, upon binding, the molecule adopts a single conformation state (thus S bound = 0), and Ω is the number of conformations in the free state, which was estimated in two different ways: i) by assigning each of the N free torsional bonds three rotational degrees of freedom (and thus Ω = 3 N ); ii) by performing a Monte-Carlo (MC) sampling with local energy minimization in the torsional space using ICM , collecting all distinct conformations within the lowest 3 kcal/mol energy, and assuming all conformers are equally probable (a similar lowlevel sampling approach was used to explore the conformational flexibility of small-molecules, Forti et al., 2012). The MC approach was considered since rotamer count is known to overestimate the number of low-energy conformations, and thus the entropy .

Quantum Mechanical Calculations
All QM calculations were performed using the QM package MOPAC2016 (Stewart, 2016) and its linear-scaling module MOZYME (Stewart, 1996), using the semi-empirical PM7 Hamiltonian (Stewart, 2013). In agreement with other authors (Sulimov et al., 2017a), we selected PM7 since it accounts for dispersion interactions, and hydrogen and halogen bonding have been taken into consideration at the paramterization stage, while it also includes several corrections to the PM6 Hamiltonian. Moreover, PM7 exhibited a very good performance on energy calculations aimed at discriminating native ligand positions in crystallographic complexes (Sulimov et al., 2017b). The solvation energy contribution in aqueous environment was calculated using the Conductor Like Screening Model (COSMO, Klamt and Schüürmann, 1993) continuum solvent model, with default atomic radii and surface tension parameters. The solventaccessible surface area was taken from the program output [cf. (Stewart, 2016) for details on how the surface is built]. Those molecules which did not complete the QM calculation were excluded when computing the enrichment.

Evaluation Metrics
The enrichment factor (EF) measures the enrichment of actual ligands in a docked hit-list given a specific percentage of the dataset (threshold). The EF is defined as the ratio between actual number ligands (hits) found at the top x% of the screened database (Hits x% ) and the number of molecules at that threshold N x% , normalized by the ratio between the total number of actual ligands within the entire dataset (Hits total ) and the total number of molecules of the latter (N total ).
EF(x) = Hits x% N x% / Hits total N total (2) Thus, the EF represents the probability of finding an actual ligand within the x% of the screened database with respect to the probability of finding an actual ligand at random. Whenever a molecule is represented within a chemical library with different states according to its protonation or chirality, each state is assigned an individual score, and the lowest score is used in the hit-list, and thus to calculate the EF. Throughout this work we report EF(1) and EF(2), since they are more representative of early enrichment. We also report receiver operating characteristics (ROC) curves for each of the studied systems, measuring the area under the curve (AUC).

THEORETICAL FRAMEWORK
The binding free energy ( G binding ) corresponding to Protein-Ligand (PL) association is expressed within the end-point molecular mechanics-quantum mechanics surface area method (MM-QMSA)  as where the difference in the first term is calculated between the bound (PL) and unbound (P, L) states, <. . . > represents the average over QM-minimized classical molecular dynamics (MD) trajectories, G QM is the QM energy including a continuum solvation term in an aqueous environment, and the second term represents the entropy change of P and L upon binding. We prefer to note the first term as a free energy, since it also includes the change in solvation free energy. Since G binding in Equation (3) is obviously too costly to be used to score and rank large chemical libraries of small-molecules in HTD, a reasonable QM docking scoring function (QMDS) can be defined as an approximation to Equation (3), namely where averages over MD trajectories have been replaced by single-point QM calculations on the docked PL structure, and the free unbound L and P structures. The L and P deformation penalty contributions due to changes in L and P conformations upon binding are expressed as where G o (X) is the energy of the isolated X in the conformation of the docked PL complex, and G(X) is the energy of X in the free unbound state. Considering Equation (5), Equation (4) can be now be written out making the deformation contributions explicit as where the "o" subscript in the first term refers to calculations using the PL, P, and L conformations from the docked complex. It should be pointed out that Equation 6 is formally identical to another formulation (Eyrilmez et al., 2019). Two types of QM docking scoring functions were defined according to the relaxation of the reference docked PL complexes: (i) QMDS 1 , with no relaxation, that is, the QM calculations are performed directly on the docked PL complex, and (ii) QMDS 2 , where docked PL complexes are relaxed through local energy minimization (see Methods). When the deformation contributions (second and third terms in Equation 6) were included, the suffix "d" is added (QMDS 1d and QMDS 2d ).

RESULTS AND DISCUSSION
Improved HTD Enrichment Using QM-Based Scoring Ten target proteins were selected based on different characteristics such as protein family, binding site properties, presence of co-factors and water molecules (within or close to the binding site), and enrichment factor at 1% calculated after docking with AutoDock Vina (Table 1). Only crystallographic and/or conserved water molecules within 4 Å of the native ligand were included. Throughout all this work, the QMDS was calculated in all its variants on PL complexes generated with ICM docking, since it is acknowledged to generate high quality protein-molecule poses (Bursulaya et al., 2003;Neves et al., 2012), as confirmed by the RMSD values of the docked native ligands in Table 2). Clearly, better enrichment is strongly coupled to scoring over correct docking poses. In this regard, the use of multiple docked conformations for each molecule, stemming from the same docking program or not, might clearly enhance the results of our QM-scoring scheme. However, we preferred to use a single pose from a single program, to keep our methodology straightforward, and to establish a clear baseline from which to start looking for improvement.
Since a target receptor protein is usually very large for QM calculations, to calculate the QMDS we used a reduced system by cutting out amino acids farther than ∼8 Å from any docked molecule (cf. the Methods section for full details on the cutout process), since a threshold of <6 Å has been reported to seriously deteriorate the results (Ehrlich et al., 2017); moreover, it should be highlighted that the smaller the threshold, the greater the impact of the continuous solvent surface replacing the cutout amino-acids. To further validate our approach, quantum mechanical docking scores QMDS considering the complete protein and its associated reduced system were calculated on CDK2 and ESR1 (Table 3). We observe that using a cutout system has no impact on the calculation. Thus, throughout this work, a reduced representation of the target protein will be used for all QM calculations.
In Table 4, we display the enrichment factors EF(1) for the 10 target systems comparing AutoDock Vina with four schemes of QM docking scoring (for HSP90a, enrichment values including and excluding the 19 macrocycle containing ligands are shown). The conformational entropy change upon ligand binding was estimated in two ways: (i) S rot , based on a term proportional to the number of N free rotatable bonds of the molecule (Ω conf = 3 N ), and (ii) S conf , by estimating Ω conf as the number of low-energy diverse conformations generated using Monte-Carlo sampling with local energy minimization (cf. Methods). We found that the use of S rot deteriorates the EF (data not shown), so S conf is used in all calculations. In QMDS 2 and QMDS 2d the reference docked PL complexes were local energy minimized using MM (see Methods). Obviously, a QM minimization would have been desirable, but this would render any QM docking scoring function useless due to the computational times involved, even for reduced systems. Moreover, in this case further caution should be exerted not to artificially deform the molecular system. As stated before, a wide range of enrichment factors calculated from docking with AD Vina was taken into account for selecting the target proteins for this benchmark. It can be readily seen from Table 4 that using any variant of QM docking scoring has an impressive improvement over AD Vina, especially in those cases with low AD Vina EF. This happens even in the simplest case of QMDS 1 , where no relaxation is performed on the PL complexes.
It is clear that PL relaxation, even using a MM-based approach, has on average a positive effect for calculating the QM docking score. Moreover, in those cases where the EF(1) slightly decreases (KITH, FABP4), the EF(2) is conserved. Focusing in the analysis of QMDS 2 and QMDS 2d , inclusion of the deformation contribution (second and third term in Equation 6) slightly deteriorates the results in ESR1, FA7, KITH, and PA2GA. However, in all but ESR1, EF(2) improves after inclusion of the deformation term (as it also happens in the other cases where EF(1) increases or is constant, CDK2, COX1, NRAM, HSP90a, HXK4, and FABP4). Considering that the effect on EF(1) is in no way dramatic, and that EF(2) (which also refers to early enrichment), improves except in one case, we state that the deformation terms are necessary to obtain better enrichment factors, though this should obviously be validated in a larger-scale benchmark. We hypothesize that this slight deterioration might be related to a small noise introduced upon energy minimization, which is canceled out in the QMDS 2 case. In the special case of HSP90a, the consideration of 19 macrocycle containing molecules has a negative effect in the EF calculation. We hypothesize that the strong performance of QM-scoring is due to a better representation of intra-and intermolecular interactions, though of course further validation and benchmarking is still needed to confirm this.
In Figure 1, the ROC plots of QMDS 2 and AD Vina for the 10 systems are shown, including the corresponding AUC values. Analysis of the curves confirm what has been noted above based on EF, exhibiting the QM-score excellent results.  Interestingly, in ESR1 both scoring methods show basically the same AUC, which is in conflict with the large difference in EF values reported in Table 4. To clarify this issue, in Figure 2 we show the enrichment plot associated to ESR1. It can be seen that AD out-performs QMDS 2 after 30% of the screened database, a region of no importance for drug discovery; for early enrichment, the enrichment plot in Figure 2 confirms the trend observed in Table 4 that QMDS 2 is remarkable superior in the initial part of the ranking. A similar behavior is observed for FA7 (cf .  Table 4 and Figures 1, 2). In the case of COX1, while the AUC of the QM-score is slightly less than AD Vina, the enrichment plot in Figure 2 shows that for early enrichment, QM-scoring out-performs AD Vina.
While our QM-score appears to be a very promising for HTD, and QM calculations are in principle more accurate than classical ones to describe molecular interactions, there are still a number of approximations which prevent the direct use of QMDS as a measure of actual absolute binding free energy. We mention three, among many: (i) QM local energy minimization was not performed (for computational efficiency, as said above); (ii) Vibrational entropies were not included; (iii) PM7 has not been parameterized to reproduce binding free energies. Our QM calculations were in the order of −70 kcal/mol, in agreement with recent binding enthalpy calculations on protein-ligand complexes using a PM7+COSMO approach (Sulimov et al., 2019a), where in spite of the difference between experimental and calculated absolute binding enthalpies, very good correlation with experimental values was obtained. It should be added that it is also well-known that traditional scoring functions correlate poorly with binding energy (cf. Enyedy and Egan, 2008, among others). Moreover, among traditional scoring functions there is no uniform scale: While AutoDock and Glide Halgren et al., 2004) are roughly in the range of −10 kcal/mol and higher, others are around −60 kcal/mol. Moreover, even end-point methods such as MM/PBSA or MM/GBSA exhibited calculated binding free energies in the order of −60 kcal/mol, or even lower when changes in vibrational entropy are not included (Zhong and Carlson, 2005), and even when including those terms (Woo and Roux, 2005;. Thus, we stress that QMDS should be considered a score, not a measure of absolute binding energy. It is aimed for relative binding energy estimation, and thus for compound ranking. On average, the computing time of this QM docking score on a single core is ∼6-8 minutes (depending on the size of the system, and on whether the deformation energy term is considered), around an order of magnitude slower than a MM-based DS.

CONCLUSIONS AND PERSPECTIVES
Docking programs have been so far based on molecular mechanics force-fields. However, a better description of protein-ligand interactions could be achieved, in principle, with quantum mechanical methods, which are theoretically exact, capture the underlying physics of the molecular system, and account for all contributions to the energy, including those effects usually missing in force-fields, such as electronic polarization, covalent-bond formation, and charge transfer. Moreover, a quantum mechanical formulation is generally valid across the chemical space, thus avoiding the force-field parameterizations.
We present a new QM-based high-throughput docking scoring function, which has been evaluated on 10 protein systems belonging to different protein families, displaying diverse binding site properties, and covering a wide range of enrichment factors computed with a traditional docking program. As shown in Table 4, even the simplest QM docking scoring function (where no relaxation is performed on the reference docked protein-small-molecule complex) shows excellent results in terms of enrichment (screening power). In fact, the improvement over AutoDock Vina on all systems is remarkable, especially in those cases with very low AD Vina enrichment. Upon complex relaxation, the improvement is even larger, regardless of whether the protein and ligand deformation terms are included or not.
We highlight that our main aim is to develop and validate a simple, straightforward approach for QM docking scoring, from which further developments can be built. Clearly, to further improve this methodology, several aspects should be analyzed: (i) a wider and extensive benchmark on many more target systems; (ii) comparison with other MM-based standard docking scoring functions; (iii) evaluation of other QM formalisms, continuum solvent models and their associated parameters (atomic radii and surface tension parameters); (iv) structural relaxation strategies; (v) use of single or multiple poses for scoring; (vi) the vibrational entropy changes upon binding. All of these considerations are important. They are currently being investigated and will be published in due course. Considering the outstanding improvements to our methods, we highlight that the QMDS should be used as a score and not an estimation to the absolute binding energy.
In terms of CPU time, our QM docking scoring function is approximately 10 times slower than MM-based standard scores on a single core. In spite of this, our impressive results on a set of 10 different protein targets highlight the huge potential of QM-based scoring. Moreover, considering future developments in QM theory, algorithms and computer hardware, it can be hypothesized that semi-empirical methods (or low-cost DFT methods) will replace FF over the following years (Grimme and Schreiner, 2018). We thus believe it is fully justified and of the utmost importance to develop the next generation of QM-based scoring functions for HTD toward highly accurate methods for the identification and optimization of small-molecule modulators of pharmaceutically relevant targets.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
CC conceived and designed the research, performed simulations, analysis, interpretation, and wrote the paper. MA performed quantum mechanical simulations and enrichment calculations and contributed with interpretation and analysis. All authors approved the manuscript for publication.

FUNDING
This work was supported by the National Agency for the Promotion of Science and Technology (ANPCyT) (PICT-2014-3599 andPICT-2017-3767).

ACKNOWLEDGMENTS
The Authors thank Prof. F. Javier Luque for helpful discussions. CNC thanks Molsoft LLC (San Diego, CA) for providing an academic license for the ICM program. The authors thank the National System of High Performance Computing (Sistemas Nacionales de Computación de Alto Rendimiento, SNCAD), the Centro de Computación de Alto Rendimiento (Computational Centre of High Performance Computing, CeCAR), and the Centro de Cálculo de Alto Desempeño (Universidad Nacional de Córdoba) for granting the use of their computational resources.