Inconsistencies in ab initio evaluations of non-additive contributions of DNA stacking energies

We evaluated the non-additive contributions of the inter-molecular interactions in B-DNA stacking by using diffusion Monte Carlo methods with fixed node approximations (FNDMC). For some base-pair steps, we found that their non-additive contributions evaluated by FNDMC significantly differ from those by any other ab initio methods, while there are no remarkable findings on their stacking energies themselves. The apparently unexpected results of nonadditivity raise issues in both FNDMC and correlated wavefunction methods. For the latter, it can be partly attributed to the imperfect complete basis set (CBS) correction scheme due to the limitation of the computational costs. On the other hand, the striking contrast between the stacking and non-additivity ∗Ryo Maezono Email address: rmaezono@mac.com (JAIST, Asahidai 1-1, Nomi, Ishikawa 923-1292, Japan) Preprint submitted to Journal of LTEX Templates October 4, 2019 ar X iv :1 80 7. 04 16 8v 3 [ ph ys ic s. ch em -p h] 3 O ct 2 01 9 behaviors was found in FNDMC. This might imply that the error cancellations of the fixed node biases in FNDMC work well for the stacking energies, while not for the non-additivity contributions involving charge transfers caused by hydrogen bonds bridging Watson-Crick base pairs.


Introduction
The non-additivity in interactions is obviously expected in inter-molecular bindings due to the induced polarizations caused by quantum fluctuations, such as vdW (van der Waals) forces. Since the binding itself has proven a great challenge to ab initio methods in terms of its description and reproduction, the further issue of non-additivity has prevented it from becoming a topic of major interest, and thus it has not been well analyzed to date. Most of the current implementations of 'molecular force fields' assume the superposition of two-body forces, applied to a vast number of simulations of self-organizations by biomolecules [1,2,3,4,5,6,7,8]. This assumption has been verified for B-DNA in a previous work [9], reporting stacking energies predicted by highly accurate quantum chemistry methods are in good agreement with those by the empirical force fields such as AMBER [10].
Recent progress in accurate computational methods, especially diffusion Monte Carlo (DMC), has enabled one to apply them to much larger systems [11,12,13,14,15,16,17,18,19,20,21,22,23]. FNDMC (fixed node DMC) [24] is a widely used implementation of DMC, where an imaginary-time dependent many-body wave function of a system evolves and gets closest to the exact (stationary) one by modifying its amplitude under the assumption that its nodal surfaces are fixed to be initial ones. The method is capable of giving an exact answer provided that the fixed node is same as that of the exact solution.
The fixed nodes are usually generated by DFT (density functional theory) or MO (molecular orbital) methods, which are not generally expected to be exact.
Hence FNDMC inherently includes a systematic bias due to the assumption on the fixed nodes (fixed node bias). Nevertheless, it has well reproduced binding natures of weakly-bound systems. [13,18,19,17] For such weakly-bound systems, it has been clarified that the non-additivity is much greater than was originally expected. [25,26] Supposing non-additive contributions were positive definite, their effects on molecular stability would get minimized in order for a system to be stable because non-additivity would evidently increase as a molecular size grows. If so, only a minor correction would be required e.g.
for C 6 (the coefficient of 1/R 6 decaying interactions) without any qualitative impact. [22] This has been demonstrated in a previous study of B-DNA stacking [9], indicating a pairwise approximation to molecular potential holds well.
Our main result in the present study cast doubt on this approximation, as shown in Fig. 4 where DMC is applied to evaluate the stacking energies of B-DNA base pairs [27,19,28,9,29,30,31] and their non-additive contributions.
While conventional methods including CCSD(T) and most of DFT predict tiny (∼ several kcal/mol), positive definite non-additive contributions, the DMC does much larger non-additive contributions with not only positive, but also negative signs, depending on base pairs; in some cases their magnitudes are as large as the corresponding stacking energies themselves (∼ 10 kcal/mol). While non-additivity in the conventional methods just slightly unstabilizes all the B-DNA base pairs, that in FNDMC can either reduce (positive sign) or enhance (negative sign) the stacking interactions depending on the pairs. The negative non-additive contributions can be deduced from a simple model analysis based on the London theory [32] as follows. Fig. 1 (b) shows the schematic geometry of Watson-Crick base pairs in B-DNA specified as 'VW:XY' in the notation convention. Base fragment pairs (W,V) and (X,Y) are respectively located at each of strands (framed by each of red-colored boxes) to form a complete four-body system. In the London theory, dispersion interactions between the upper and lower layers can be approximated as ε ∼ α (upper) · α (lower) , where α (upper) and α (lower) are the dipole polarizabilities of each layer. Note that the polarizability is roughly additive [33] because it scales as the molecular weight. This makes a rough estimate of the stacking energy for the entire (four-body) system as ε (4) = 2 × 2 · ε (2) < 0, where ε (2) denotes the stacking energy for a partial (two-body) system. (Of course, there are other dependencies of ε such as on ionic energies and geometry etc., which do not have such a great effect -our data given in the Supporting Information have demonstrated that the simplified discussion holds.) The estimate then gives the nonadditivity as ∆ε (4) = ε (4) − 4 × ε (2) = 0. This would be true only for the limit, l → 0 [(a /a) → 1]. In practical cases [(a /a) > 1], however, we can ignore "cross" interactions [between W and Y and between X and V (see Fig. 1)] due to (1/a 6 ) (1/a 6 ) ((a /a) = 1.5∼1.9). We then arrive at the negative non-additivity: ∆ε (4) Upon applying the simplified model analysis, the London theory produced a negative sign of the non-additivity contribution, implying that the stacking energy is enhanced by the non-additivity.
It is worth noting that the negative non-additivity found in FNDMC and the above-mentioned model analysis is quite different from the results obtained using other ab initio methods. In this paper, we provide several possible theories explaining why discrepancies arise between the FNDMC and the other predictions. We first investigated that 'negative values are predicted only by DMC'.
Once we verified this issue, a further doubt occurs regarding the sign alternation.
We also provide plausible explanations for the sign of non-additivity.
Our central findings in the present study are summarized as follows: Our evaluation of binding energies come up to common expectations for the methodologies (Fig. 3). But this appears not the case for non-additivity at first glance (2) "CCSD(T)" gives almost the same results as B3LYP and Hartree-Fock (HF), which hardly happen when evaluating the binding energies. These apparently strange conclusions are likely to mislead readers, so we shall out-line the contents below. We claim that the conclusions are not due to careless choices of our computational specifications, but more fundamental points -practical approximations adopted in DMC and CCSD(T), each of which works well for evaluating binding energy, but not for non-additivity. In the former case, we can see many previous works reporting the cancellation of the fixed node biases between a whole system and its constituent molecules works well for evaluating binding energies [18,19,17]. On the other hand, such a cancellation has not been investigated yet. In the latter case, we note the fact that "CCSD(T)" applied to B-DNA systems is actually "CCSD(T) with CBS at the MP2 level" [9].
We finally conclude this practical approximation can be attributed to the reason why it gives the same trends in non-additivity as B3LYP that is not believed to be capable of reproducing vdW interactions [18].

System and methods
Our target systems are ten unique B-DNA (Watson-Crick) base-pair steps shown in Fig. 1. Their molecular geometries were obtained from a pioneering work due toŠponer et al. [9], and they are fixed throughout our simulations.
Note that, e.g., 'AT:AT' and 'TA:TA' seem identical schematically because of their mirror image relation. Actually, they have different geometries in B-DNA sequences, and hence they are not identical. Stacking energies, ε (4) and ε (2) are evaluated using ab initio methods. In the present study we adopted the same non-additive contribution ∆ε (4) as that defined byŠponer et al. [9]: VW +ε VX +ε Note that this definition does not include the contributions of the hydrogen bonding between bases in a Watson-Crick pair within each step: V:Y and W:X.
The dispersion interaction -a main contribution to stacking -can be reproduced by going beyond MP2 (Møller-Plesset) level treatment of electron correlations [26,28,34,35]. We hence applied FNDMC (diffusion Monte Carlo with the fixed-node approximation) [36] to compute the ε (4/2) values in Eq. (1) !""""#"""""$ !""""#"""""$ %&' %(' using the CASINO code [37]. Computational details in our FNDMC were the same as those in our previous studies [18,20]. Furthermore, they met a protocol developed for non-covalent systems by Dubecky et al. [38]: We adopted Burukatzki-Filippi-Dolg pseudopotentials (BFD-PPs) [39] and applied the Tmove scheme [40] to their evaluation with the locality approximation [41]; fixednode guiding functions were constructed from DFT-B3LYP/VTZ orbitals with Gaussian09 [42]; one-and two-body Jastrow functions were considered where their parameters were optimized by variance minimization procedure [43,44]; a time step is set to be δτ = 0.005, which is small enough to avoid the time step error bias. [45] A number of previous works have demonstrated that FNDMC with similar conditions are computationally feasible for non-covalent systems as large as the present B-DNA systems and successfully reproduce their non-covalent interactions comparable to CCSD(T)/CBS (coupled-cluster approach at the singles and doubles level augmented with perturbative triples using complete basis set limit). [11,12,13,14,15,16,17,18,19,20,21,22,23] In particular, the B-DNA stacking energies evaluated by FNDMC were calibrated in detail in our previous works. [18,20] There is no doubt that CCSD(T)/CBS is a state-of-the-art or "gold standard" quantum chemistry method -an established protocol of capturing dispersion interactions in non-covalent systems. [46] Its applicability is, however, quite limited to small systems. Actually, CCSD(T)/CBS is not applicable to our target systems of B-DNA base-pair steps, even using the best supercomputer. The practically best possible solution [9] was to apply "CCSD(T)/CBS" to all pairwise stacking and add a many-body correction at MP2/VDZ level to the pairwise sum. Here we note that "CCSD(T)/CBS" for the pairs is not a true one, but an approximation such that MP2/CBS is combined with a energy difference between CCSD(T) and MP2 obtained using a small basis set.
Hereafter we refer this sort of approximation to CCSD(T)/CBS(MP2). Very recently, Kraus et al. [47] have attempted to avoid the approximation of stacking interaction such as the sum of four base-base stacking energies, but their level of theory has not reached CCSD(T)/CBS. Parker et al. [31] stated that such an approximate estimate can be used for reference, but not as a conclusive standard value. In the above-mentioned context, a true non-additivity at CCSD(T)/CBS level of theory remains unknown.
In addition to the 'correlation-level' non-additivity, we also investigated 'SCF-level' non-additivity by performing Hatree-Fock (HF) and DFT with various choice of exchange-correlation (XC) functionals. Our HF and DFT simulations were carried out using Gaussian09 [42] with the same basis set and pseudo potential as FNDMC. It is well known that some of our choices of XC functionals here (e.g., LDA and B3LYP) are inappropriate for describing the interactions in the present system, but we intentionally adopted them for comparison between correlation-and SCF-level non-additive contributions. In addition to the conventional XC functionals, we adopted recently developed XC functionals such as ωB97X, ωB97M-V [48], B3LYP-D3 [49], and CAM-B3LYP-D3 [50] in order to investigate not only dispersion effects, but also hydrogen bonding.

Results
The present study adopts CCSD(T)/CBS[MP2] results due toŠponer et al [9] as reference, though they are not "true" CCSD(T)/CBS as described in the previous section "System and methods". Here we show not only nonadditive contributions ∆ε (4) , but also stacking energies ε (4) themselves similar to conventional researches. We will see that FNDMC gives rise to a striking dependence of ∆ε (4) on base-pair steps, compared to the other ab initio methods including CCSD(T)/CBS[MP2] as well as DFT, while there is not any novel findings in their ε (4) evaluations.

Stacking energies ε (4)
We shall start with the four-body stacking energies (ε (4) ) of the B-DNA basepair steps. as reference. Within the wavefunction methods ( Fig. 3 (a)), we can see that HF fail to describe the stacking properly because it does not include electron correlations at all by nature. It is obvious that an appropriate description of stacking requires theory more than MP2 level [34,35]. Overall, the correlated methods agree with each other. Looking closely at the trend, CCSD(T)/CBS [MP2] was found to overbind compared with the other approaches, LMP2, SAPT and FNDMC. This can be attributed to two practical approximations adopted in the CCSD(T) reference [9]: Since the B-DNA system is too large to be computed by CCSD(T) with larger basis sets (cc-pVTZ and cc-pVQZ), the true CCSD(T)/CBS is not available in any literature so far. Instead, two approximations were applied to estimate "CCSD(T)/CBS". [9] Firstly, all the pairwise base-base terms (ε (2) in Eq. (1)), were computed by MP2/CBS plus an energy difference between CCSD(T) and MP2 with a common small basis set (6-31G*(0.25)). Secondly, the non-additive contribution (∆ε (4) ) was evalu- Step overbinds, similar to MP2. Note that both SAPT and DF-LMP2 are known to correct the overbinding trend in MP2 [29]. Accordingly, 'the overbinding ' can be also supported by the fact that the magnitudes of stacking energies in SAPT and DF-LMP2 are smaller than those in , as can be seen in Fig. 3 (a). It was found that FNDMC deviates from the other correlated methods depending on base-base pair steps (especially TA:TA and AG:CT), which is closely related to a significant difference in non-additive contributions between the methodologies, as described later.
As for the DFT methods, we see that only B3LYP cannot properly describe the stacking, which is consistent with its well-known deficiency in describing dispersion interactions [51]. which has been also demonstrated for other non-covalent systems [19]. Note that the stacking described by LDA is known to be artificial. [13,14,18,19,20,22,52,53] 3.2. Non-additivity ∆ε (4) Non-additive contributions ∆ε (4) evaluated from various methods are shown in Fig. 4. We found FNDMC giving a wiggling behavior of ∆ε (4) , compared to all the other methods. This remarkable sign alternation found in FNDMC is a central issue in the present study.
We first remark that the non-additive contributions do not necessarily arise from the electron correlations but always appear as non-linear processes inside a many-body system. In terms of the perturbation expansion based on SAPT, an interaction energy is decomposed into physically meaningful components: electrostatic, induction, exchange, and dispersion terms. [34,35] Even at the Hartree-Fock (HF) level of theory, the 'exchange' and 'induction' parts of the non-additivity occur [34,35], which we refer to 'SCF-level non-additivity' hereafter. The behavior of HF in Fig. 4 can be an appropriate reference to the 'SCFlevel non-additive' contribution because HF is incapable of describing the dis-  A comparison between B3LYP and B3LYP-D3 gives the most suggestive insight into the non-additivity within the framework of DFT. While the empirical dispersion correction D3 significantly improves the stacking itself ( Fig. 3 (b)), it hardly modifies its non-additivity from B3LYP ( Fig. 4(b)). This can be attributed to the fact that dispersion corrections based on D3 or the likes of vdW-XC are additionally made on the original DFT/SCF energies and thus, never deform their wavefunctions. From the viewpoint of many-body theory, the deformation of wavefunctions is the origin of dispersion interactions and hence essential to the non-additivity.

Sign alternation of ∆ε (4) in FNDMC
Except for B3LYP (-D3)  In summary, the dispersion contributions to non-additivity have the potential for changing the sign of non-additivity, depending on base-pair steps.
If we accept the negative non-additivity in FNDMC in accordance with the London theory, then another doubt arises about why FNDMC also gives more positive non-additivity, depending on the base-pair steps. Comparing FNDMC with HF, FNDMC was found to give almost the same ∆ε (4) as HF (within errobar) for GG:CC, GC:GC, and CG:CG (Fig. 4). In contrast, the more positive ∆ε (4) deviating from the SCF-level non-additivity appears in TA:TA, GA:TC, and AG:CT, which will be then investigated from the viewpoint of stacking energies. Fig. 5 shows one four-body and four two-body stacking energies given in Eq. (1) for ten unique B-DNA base-pair steps. It is evident that positive nonadditivity correlates with weaker four-body stacking (red bar), which is also noted in Fig. 3 that the base-pair steps with the positive non-additivity exhibit the weaker stacking described by FNDMC than by the other ab initio methods. This factor was not taken into account in our simple London model analysis. To estimate the contributions quantitatively, we first evaluated the Mulliken charge that appeared in the bridging location, as shown in Fig. 6 (c). Based on the charge, we then evaluated the Madelung repulsion interaction, getting +5∼10 kcal/mol per bond. Since the typical range of the energy gain due to the hydrogen bonding is known to be less than 5∼6kcal/mol [56], it is likely to result in a positive contribution to non-additivity, thus making the stacking weaker.

Discussion
The present results are unexpected against a common consensus among the methodologies, and hence likely to mislead readers. Here we shall clarify what points are to be noted in order to avoid misleading. The first point is that the present results never immediately lead to such a simple controversy as 'DMC vs CCSD(T)'. The second point is that the unexpected results are obtained only at the non-additivity level (Fig. 4), but not at the stacking energy level (Fig. 3).
As for the first point, our results in Fig. 4 are apparently regarded as 'DMC vs [CCSD(T) together with all the other DFT/SCF methods]'. This might be followed by such a misleading discussion as 'which method is reliable?', but it is actually premature to draw a conclusion about this question. We should clarify limitations on approximations adopted in the present FNDMC and CCSD(T) simulations.
We start with CCSD(T) together with HF/DFT. It was seen that HF and B3LYP failed to describe the B-DNA stacking interactions ( Fig. 3 (a) and (b)), but their non-additive contributions were almost same as "CCSD(T)" (Fig. 4 (a) and (b)). This must be unacceptable because HF and B3LYP do not include dispersion interactions at all. [34,35,49] (4) as the other DFT methods. In particular, the dispersion correction based on B3LYP-D3 significantly improves the staking itself, but seldom changes the non-additive contributions from B3LYP. This means that reproducibility of ε (4) hardly affect that of ∆ε (4) at SCF level. MP2 can reproduce a main part of the dispersion interaction (or electron correlation) by the second-order perturbation energy correction, but its wavefunction still remains the Hartree-Fock one. This is analogous with the relation B3LYP and B3LYP-D3. Thus, we may conclude that the dispersion correction accompanied with wavefunction deformation is essential for reproducing the dispersion-level non-additivity.
The overall coincidence between "CCSD(T)" and SCF (HF/DFT) implies that (1) the present "CCSD(T)" method never describe the dispersion-level non-additivity and thus, losing a large part of true non-additivity at CCSD(T) level of theory, or (2) a true dispersion-level non-additivity is essentially tiny and thus, well described by the present "CCSD(T)" method. As has been (as well as the true CBS) to the B-DNA base-pair steps, but they were too large to compute. Instead, we dealt with a neon tetramer as a simple/model system in which each DNA base is replaced by a Ne atom. We evaluated ∆ε (4) values of the Ne tetramer at several distances between the two dimers using  Figure 7: The non-additivity ∆ε (4) of neon tetramer at several distances between the constituent dimers (described as "Interlayer distance" in the horizontal axis). With a fixed "Interlayer distance", all the Ne atoms located on a plane form a rectangle, where in each dimer its interatomic length is fixed to be 2.925Å. All the CCSD(T) and DFT calculations were performed using Gaussian09 [42].
Next we move on to FNDMC. Its wiggling dependence of ∆ε (4) on the basepair steps appearing in Fig. 4 arouses suspicion whether FNDMC really reproduces the dispersion-level non-additivity appropriately. Here we shall deliberate on two possibilities of causing faults in FNDMC: quality of trial wavefunctions obtained and reliability of the fixed-node approximation adopted in the present study.
We first note that such a wiggling dependence appears only in Fig. 4, but not in Fig. 3; both the results were obtained by the wavefunctions that were optimized at the same level of theory. In the present study we did not optimize the Slater part the trial wavefunctions, but the Jastrow part only. In FNDMC, the latter changes the statistical error bar only, while the former -related to the fixed-node approximation -changes the final total energy value, unlike VMC (variational Monte Carlo). [24] It is well known that the same performance on the Jastrow optimization leads to the same magnitude of error bars for similar system sizes if all the other computational details are assumed to be common to the systems. It is obvious from Figs. 3 and 4 that this is actually valid for the present B-DNA systems. We also insist that our choice of computational details on FNDMC -basis sets, time step, t-move scheme, as well as Jastrow function -is equivalent to a protocol established in previous studies on noncovalent systems due to Dubecký et al. [17]. The reasonable behavior of FNDMC stacking energies in Fig. 3 asserts that our choice is valid for the present B-DNA base-pair steps.
The fixed-node approximation is the most notorious as the cause of errors in FNDMC. [36] Previous studies on non-covalent systems including B-DNA, however, demonstrated that FNDMC works well for evaluating their complexation energies in general. [17] It is to be noted in the B-DNA stacking that this is valid for the stacking energies, but unknown for the non-additivity. The success in the FNDMC stacking energies relies on the error cancellation of the fixed-node approximations between the whole non-covalent system and its constituent subsystems. This implies that the formation of non-covalent/vdW bonding does not give rise to a significant difference in nodal surface structures between the whole and the sub-systems (tetramer-dimer/dimer-monomer), leading to an accurate complexation energy. The success in the non-additive contribution requires that two error cancellations of the fixed-node approximations simultaneously occur for ε (4) and ε (2) . In the case of non-additivity, however, it is possible that the cancellation in ε (4) would not occur properly. Its possible factor can be attributed to a horizontal bridging between Watson-Crick bases due to hydrogen bonding. The formation of hydrogen bonding accompanied by the charge transfer could deform the fixed-node surface structure of the tetramer more significantly than that of vdW bonding. If this were true, the fixed-node error could not be canceled out more remarkably in ε (4) than in ε (2) . The less cancellation could arouse the suspicion that the wiggling dependence of FNDMC in Fig. 4 is incorrect due to the fixed-node error related to the hydrogen bonding.
In order to prove the conjecture (or anti-conjecture) about the fixed-node errors caused by the hydrogen bonding, the most straightforward way would be to evaluate the nodal surface dependence of the non-additivity. While the stacking energy has been demonstrated to be insensitive to the dependence [18,20], one would suspect that it is not the case for the non-additivity. Suppose it were true, the non-additivity evaluated with a different trial node would be different from the present FNDMC one in Fig. 4. Although we plan to address this issue in our future work, we should mention that such a calculation involves a heavy computational resource, which costs 1.2×10 6 core-hour [1.2×10 5 (corehour) × 10 pairs]. In addition to a single reference trial node, we could employ recently developed trial nodes such selected configuration interaction [57] and multipfaffian [58] wavefunctions, as well as a simple multi-reference one [59,23].
Beside their feasibility in terms of computation costs, the more sophisticated trial nodes would shed light on the nodal surface dependence of non-additivity in FNDMC, i.e., we could verify whether or not the non-additivity is sensitive to the nodal surface unlike the stacking energies [18,17].
Although we do not investigate the nodal surface dependence of ∆ε (4)   Watson-Crick base giving rise to a large fixed-node error, and hence the wiggling dependence of ∆ε (4) in FNDMC, shown in Fig. 4, might be false due to the less error cancellation. in the ∆ε (4)   we could not obtain numerically/statistically reliable FNDMC results because the magnitude of ∆ε (4) itself is an order of/less than the sub-chemical accuracy (0.1 kcal/mol) and hence the corresponding error bar is required to be an order of/less than 0.01 kcal/mol. To attain such an error bar, a vast number of statistical samplings must be accumulated even for the smaller system considered here. This is another serious drawback to be noted for FNDMC.

Conclusion
We applied a diffusion Monte Carlo method with the fixed-node approximation (FNDMC) to evaluate the non-additive contribution to B-DNA stacking.
We found that the FNDMC values of non-additivity alter their sign (i.e. they increase or decrease their stacking interactions) depending on the base-pair steps, which is contrary to all the other ab initio methods. On the other hand, no significant difference between the methodologies was observed for four-/two-body stacking energies, each of which are used to evaluate the non-additivity. To elucidate this contrast between the stacking and non-additivity, we made two plausible discussions about limitation on practical approximations involved in CCSD(T) and FNDMC: (1) The reason why the unexpected coincidence between CCSD(T)/CBS[MP2] and HF/B3LYP occurs only at the non-additivity level can be attributed to the imperfect capability for MP2 to reproduce the electron correlation specific to the four-body system. The lack of the correlation never describes the correlation/dispersion-level non-additivity properly. In other words, CCSD(T)/CBS[MP2] mostly describes the SCF-level non-additivity only.
(2) FNDMC demonstrates a wiggling dependence of the non-additivity. While the SCF-level non-additivity is mostly positive, the non-additive contributions described by FNDMC are both positive and negative signs. The negative sign is found to be reasonable, which might be supported by a simple model analysis based on the London theory. It would, however, be premature to draw a conclusion that the FNDMC non-additivity reveals the truth. This is because the Watson-Crick base-pair involves the charge transfer caused by the hydrogen bonding, but we could not verify if the error cancellations of the fixed-node errors were successful for the hydrogen bonding, as in the case of complexation energies. In the near future, we will investigate the fixed-node dependence of the non-additivity by considering more reliable trial nodal surfaces, albeit with huge computational costs.

Acknowledgments
The computation in this work has been performed using the facilities of the