{'fraction_non_alphanumeric': 0.0254411161263849, 'fraction_numerical': 0.007509232663110382, 'mean_word_length': 4.363336267605634, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Simple and Effective Solution to the Constrained QM/MM Simulations

It is a promising extension of the quantum mechanical/molecular mechanical (QM/MM) approach to incorporate the solvent molecules surrounding the QM solute into the QM region to ensure the adequate description of the electronic polarization of the solute. However, the solvent molecules in the QM region inevitably diffuse into the MM bulk during the QM/MM simulation. In this article we developed a simple and efficient method, referred to as boundary constraint with correction (BCC), to prevent the diffusion of the solvent water molecules by means of a constraint po- tential. The point of the BCC method is to compensate the error in a statistical property due to the bias potential by adding a correction term obtained through a set of QM/MM simulations. The BCC method is designed so that the effect of the bias potential completely vanishes when the QM solvent is identical with the MM solvent. Furthermore, the desirable conditions, that is, the continuities of energy and force and the conservations of energy and momentum, are fulfilled in principle. We applied the QM/MM-BCC method to a hydronium ion in aqueous solution to construct the radial distribution function(RDF) of the solvent around the solute. It was demonstrated that the correction term fairly compensated the error and led the RDF in good agreement with the result given by an ab initio molecular dynamics simulation.


I. INTRODUCTION
The hybrid quantum mechanical / molecular mechanical (QM/MM) approach has been extensively utilized to study solvation processes in chemistry and physics. [1][2][3][4][5][6][7][8][9][10] Actually, it offers a versatile theoretical framework for describing the electronic state of a QM object under the influence of a surrounding MM environment. The polarization of the electron density of the QM region due to the environment is a matter of a critical importance in the QM/MM simulations. The success of the QM/MM method in condensed phase simulations is attributed mainly to the fact that the polarization of the QM solute can be reasonably realized by the Coulomb potential due to simple point charges placed on the MM molecules.
We note, however, that the charge-transfer type polarization between QM and MM objects cannot be considered in the QM/MM method at least in a theoretically proper manner.
In particular in a hydration of an ion, where the charge-transfer type polarizations play a role, the QM/MM approach is no longer adequate to evaluate the solute-solvent interaction in general. Furthermore the orbital mixings due to the intermolecular charge migration will seriously affect the exchange repulsion 11,12 between solute and solvent molecules defined usually by the size parameters in the Lennard-Jones potentials 13 . It is, thus, desirable to develop a methodology to resolve the problems associated with the QM/MM boundary to extend the frontier of the QM/MM approach.
A straightforward but promising solution to the problems is to incorporate the solvent molecules surrounding the solute of interest into the QM region. This treatment offers an adequate description of the electronic states of the solute interacting with the neighboring solvent molecules. However, such a preferred solvent configuration will collapse during a simulation due to the free diffusion of the QM solvent into the MM bulk solvent. To prevent the diffusion of the QM solvent, there have been a lot of developments that are classified into two categories referred to as 'adaptive' [14][15][16][17][18][19][20][21] and 'constrained' 22,23 QM/MM approaches.
In the former approaches, the solvent molecules can change its character adaptively from QM to MM (and vice versa) when they cross the QM/MM boundary. Most of the adaptive approaches introduce a transition(buffer) zone between QM and MM regions to ensure the continuity of the energy and forces of the solvent. In the latter, on the other hand, some constraint forces are applied to the solvent to keep the solvent molecules within a solvation shell. We note, however, that any constraint or adaptive procedure violates a principle of the statistical mechanics since the situation that some particular solvent molecules stay within a shell is obviously against the entropy effect. In other words, the mixing entropy will be inevitably lost in a constraint approach. Thus, the requirement we can impose at most on a boundary treatment is that the artifact completely vanishes when the force field of the QM solvent is identical with that of the MM solvent. This may constitute a principal condition to be fulfilled in the development of the adaptive and the constrained QM/MM methods.
Actually, the adaptive QM/MM approach will satisfy this condition since it is naturally expected that the effect of an adaptive procedure will disappear when the QM solvent is identical to the MM solvent. Rowley 23 which applies bias force to every pair of the QM and MM particles to prevent the two regions from being mixed. In practice, however, the bias forces will be applied only to the pairs with significant contributions to reduce the large computational costs.
In the present work we develop a new method of constrained QM/MM, called Boundary Constraint with Correction (BCC), which fulfills the desired conditions. That is, the continuities of energy and force, and the conservations of energy and momentum are guaranteed in the BCC method and more importantly BCC is designed so that the effect of an applied bias potential completely disappears when QM solvent is identical to MM solvent.
Furthermore the BCC method can be readily implemented provided some QM/MM code is available. As will be described in the Theory section the basic idea underlying the BCC approach is very simple. The point of the method is to make a correction compensating the error in a statistical property of interest due to the constraint force by performing a set of separate QM/MM simulations. Thus, the BCC method necessitates an additional procedure in contrast to the FIRES 22 and BEST 23 methods. However, the computational load associated with the correction is rather modest. The efficiency of our method is first assessed by conducting full MM simulations for water solutions where the water solvent is partitioned into two regions through a bias potential. Then, the radial distribution functions (RDF) 13 obtained by BCC are compared with those yielded in conventional simulations without the bias potential. The robustness of the method is also examined by varying the number of water molecules contained in each domain. We next utilize the BCC method to study the hydration of a hydronium ion in water solution where H 3 O + and water molecules within the first solvation shell are represented by a density functional theory(DFT) for electrons. The RDF for water solvent around H 3 O + is constructed through a QM/MM simulation combined with BCC.
This article is organized as follows. In the next section we show the theoretical framework of the BCC approach that can be fully expressed with a single equation. In Section III the computational details for the full MM and the QM/MM simulations are presented. In Section IV the RDFs given by the MM simulations combined with BCC are compared with those obtained by a conventional approach to discuss the efficiency and the robustness of the BCC method. The hydration structure for H 3 O + constructed using the QM/MM−BCC method is also compared with a recent result provided by a first-principles simulation. Finally, in Section V we summarize our work and make a remark on the prospect of the extension of BCC to free energy calculations.

II. THEORY AND METHODOLOGY
This section consists of two subsections. In Subsection A the details for the boundary constraint with correction (BCC) will be described. Then, the explicit form for the constraint potential employed in the present work will be given in Subsection B.

A. Boundary constraint with correction
To introduce the basic idea of the BCC approach we consider here a solution described with pairwise potentials. The discussion can also be extended to a QM/MM system straightforwardly. We suppose that the solvent consists of two kinds of molecules A and B with different solute-solvent and solvent-solvent potentials. As illustrated in Fig. 1 solvent A is being confined within a solvation shell Ω enclosing a solute. Then, the number N A of molecules A would be much smaller than N B for the bulk solvent B in general. For our present purpose it can be assumed that the force field of molecule A is reasonably close to that of B. Provided that x A collectively denotes the configuration of the solvent molecules A the interaction energy U A (x A ) can be given by where u A and υ A are, respectively, the solute-solvent and solvent-solvent interaction potentials for solvent A, and x Ai denotes the full coordinates of ith molecule. In Eq. (1) it is assumed that the position and the orientation of the solute molecule is being fixed during a simulation. We also have an equivalent equation for solvent B, thus, The interaction U AB (x A , x B ) between solvent A and B is written as As illustrated in Fig. 1 we apply a bias potential U bias to solvent molecules A to keep them staying within a solvation shell Ω. In terms of these potential energies, the statistical average bias of a physical property P under the constraint potential can be expressed by In Eq. (4) β is the inverse of the Boltzmann constant k B multiplied by temperature T , and the coordinate X denotes the set of the coordinates x A and x B , i.e. X = (x A , x B ). We are interested in how to eliminate the contribution due to the bias potential from the statistical To this end we consider the solution with 'homogeneous' solvent in which solvent A also obeys the potentials (u B , υ B ) in Eq. (2), and υ AB in Eq. (3) is identical to υ B . Then, the ensemble average P u B ,υ B bias under the bias potential is written as For the homogeneous solvent system, the average P u B ,υ B without the bias potential is also well defined and merely given by The quantity of Eq. (6) subtracted by Eq. (5) can be regarded as the correction to the physical property P to compensate the effect due to the boundary constraint exerted on molecules A. We, thus, define the correction term ∆ corr P u B ,υ B as Therefore, Eq. (4) will be corrected as It is readily recognized in Eq. (8) that the effect of the bias potential completely disappears when u A = u B and υ A = υ B = υ AB . Equation (8) describes the framework of our boundary constraint approach. As a major application of our method we consider a QM/MM system where the solvent A as well as a solute molecule are described with a quantum chemical theory and the solvent B is expressed with an MM force field. In such a system, Eq. (8) can be rewritten as Of course, in Eq. (9), E QM and E MM are the energies of the QM and MM regions in the system, respectively, and E QM/MM is the interaction between the two regions.
In closing this subsection we make a remark on the BCC approach. Since the correction in BCC is yielded through the statistical averaging as shown in Eq. (7) it is not possible to apply it to evaluate some dynamical properties. We also note that the solvent molecules B are allowed to enter the solvation shell since no constraint forces are directly applied to them. However, the solvation shell Ω will be exclusively occupied by the solvent A and the effect of the contamination due to solvent B will be marginal. Anyway, we should determine the number N A beforehand for the given size of Ω. The effect of the choice of N A will be examined in Results and Discussion. As shown in Eq. (8) the BCC energy consists of three ensemble averages, each of which utilizes a trajectory yielded in a continuous molecular dynamics simulation conserving energy and momentum.

B. Potential for constraint
In this Subsection we introduce the bias potential U bias employed in the present work to constrain the solvent A within a solvation shell Ω. Actually, the choice of U bias is somewhat arbitrary as long as it well behaves. Our choice in this work is to adopt the logarithm of the Fermi function that was proposed in Ref. 23. The explicit form of U bias is given by where r is the solute-solvent distance and r c is the radius of the spherical shell Ω assuming that the center of a spherical shell Ω is placed at an interaction site of the solute. It can be easily recognized that U bias (r) is almost zero in the region of 0 ≤ r ≤ r c , while the slope of U bias (r) asymptotically approaches k B T × α for r c ≤ r. Thus, a solvent molecule feels no constraint force when it is inside the shell Ω, while it will be repelled by a slope when it exits the shell. A half harmonic potential can also be used as a bias potential U bias though not tested in the present work.
As described above the shape of the shell Ω is assumed to be spherical in the present benchmark calculations. In a practical application, however, a more complex shape will be appropriate depending on the structure of the solute molecule of interest. An instant solution to this problem is for example to make a union of the spherical shells {Ω i } (i = 1, 2, · · · ), each of which center is placed at the ith interaction site of the solute. Of course, such a non-spherical shell can be readily adopted in the BCC method. We note a multiple-sphere scheme was also proposed in the adaptive QM/MM method. 16

III. COMPUTATIONAL DETAILS
In the following we apply the BCC method to two kinds of aqueous solutions as benchmark calculations. We first consider a solution in which a water molecule is solvated by water solvent. All of the molecules in the system are represented with classical force fields. 13 The solvent consists of two kinds of molecules A and B specified with different potentials. A constrained potential is applied to solvent molecules A to keep them within a solvation shell Ω around the solute. Then, the radial distribution function (RDF) 13 for the solvent will be obtained through Eq. (8). The details of the simulations are provided in Subsection A.
Next, the hydration of a hydronium ion(H 3 O + ) is studied by a QM/MM simulation combined with the BCC method, where a hydronium ion and the neighboring 4 water molecules are involved in the QM region through a bias potential. We also compute the RDF for water molecules around the solute to make comparisons with a reference result given by an ab initio molecular dynamics (AIMD) simulation. 24 The computational details for the QM/MM-BCC simulations are given in Subsection B.

A. A water molecule embedded in a water solvent
We apply the BCC method to solutions described with fully classical force fields. In each solution system, a water solute is embedded in a mixed solvent comprised of two kinds of water molecules A and B with different interaction potentials. In applying the BCC method, solvent A is confined within a shell Ω through a constraint force whereas no bias potential is applied to solvent B as illustrated in Fig. (1). For all the water molecules in the system the force field parameters except for the point charges refer to the SPC/E model. 25 Solvent B is completely identical to SPC/E model and has a negative charge q B o = −0.8476 in the unit of the elementary charge on the oxygen site, while in solvent A it is shifted to q A o = −0.90 assuming a larger polarization. The solute molecule is assumed to be identical to solvent A and hence the charge q o on the oxygen site is taken as −0.90. Of course, the charge on a hydrogen atom is determined to ensure the neutrality of the water molecule. To assess the efficiency of the BCC approach we construct the radial distribution functions(RDF) through Eq. (8) for the solvent oxygen around the solute oxygen. The Newtonian equations of motion for the solvent molecules are numerically solved using the leap-frog algorithm 13 and the internal structure of the solvent molecule is fixed during the simulation using the quaternion. 13 The time step for the molecular dynamics(MD) is set at 1 fs. The solute water molecule is fixed at the center of a cubic simulation box with a periodic boundary condition. 13 Total number of the water molecules contained in the simulation box is 500 including the solute. The long-range intramolecular Coulomb interaction is evaluated using the Ewald method. 13,26 The thermodynamic condition is set at temperature T = 300 K and density ρ = 1.0 g/cm 3 . The statistical average for each RDF is yielded from a trajectory of 2 ns MD simulation.
We conduct the BCC simulations employing the bias potential defined by Eq. (10) The Lennard-Jones parameters 13 for SPC/E is assigned to the oxygen atom of the hydronium ion. The RDF for the solvent oxygen around the solute oxygen is constructed to make comparisons with that given by the first-principles approach. 24 The molecular geometry of the solute is optimized by the KS-DFT with the B3LYP exchange correlation functional [29][30][31] and the aug-cc-pVTZ basis set. 32 The optimization is performed by conducting Gaussian 09 program package. 33

IV. RESULTS AND DISCUSSIONS
This section is also partitioned into two subsections in parallel to the organization of the Section 'Computational Details'. Subsection A is devoted to show the results of the BCC simulations for a water molecule immersed in a water solvent, and in Subsection B the results and discussions are provided for the QM/MM-BCC simulations applied to a hydronium ion in water solution.  In our constraint method solvent B is allowed to enter the inside of the shell Ω. Hence, the ratio of the occupancy of the shell by solvent B is also of our concern. To this end we computed the individual RDFs for solvents A and B. The results are depicted in Fig. 2(b) which also shows the bias potential of Eq. (10) as a function of the O−O distance r. It is seen in the figure that the inner region of the shell is partly occupied by solvent B. Explicitly, the ratio of the occupancy of solvent A to B is estimated as ∼ 5 at the first peak of the RDF.

A. A water molecule embedded in a water solvent
Thus, it was revealed that solvent B also enters the inner region of the shell although solvent A dominates the region. This feature constitutes a difference of our approach from other constraint methods such as FIRES 22 and BEST 23 .
It is of fundamental interest how the bias potential affects the orientational structure of the solvent molecules around the solute water. To quantify the tetrahedral order we introduce the orientational order parameter q defined in Ref. 42, thus,  We also applied the BCC method to a system with larger shell Ω(r c = 7.5Å), for which N A was set at 58. The RDFs relevant to the BCC calculation is shown in Fig. 3(a) and then, the number N av was taken as N A . We note, however, that the instantaneous number N in of the molecules inside the shell substantially fluctuates during the simulation.
Actually, the minimum and the maximum numbers of water molecules in the shell were 48 and 69 in our simulation, respectively. Therefore, confining N av solvent molecules constantly within the shell throughout the BCC simulation will give rise to some artifacts. It is, thus, suggested that the number N A for a given size of the shell should be chosen so that it allows the fluctuation of the number of the water molecules contained in the shell during a simulation. In other words, N A should be appropriately smaller than N av considering the deviation of N in . The choice of N A = 50 is, thus, found to be suitable for the wall position of r c = 7.5Å. The calculation for N A = 48 was also performed and the result is plotted in Fig.  3(b), which shows little difference with the plot for N A = 50 as expected. We speculate that the naive choice of N A = N av for the smaller shell (i.e. r c = 3.5Å) gave a rather successful result since N in does not fluctuate largely. Anyway, the proper number N A for a shell Ω can be readily determined through a preliminary calculation.
In closing this subsection we make a brief remark on the correction term of Eq. (7). As demonstrated in the above simulations, the correction works adequately when the potential of A is reasonably close to that of B. However, it is possible that a flaw will emerge in the RDF when the force fields for A and B are largely different from each other. Actually, the artifact of the bias potential remains in principle when solvent A is not identical to B in the BCC approach as in the other constraint approaches. In the following subsection we also test the BCC correction by applying it to an actual QM/MM system involving an ionic molecule as a solute of interest.

B. QM/MM-BCC simulation for hydronium ion
We combine the QM/MM simulation with the BCC approach (QM/MM-BCC) to study the hydration of a hydronium ion. The H 3 O + and 4 water molecules confined within a shell Ω (r c = 3.5Å) were described by the KS-DFT method. The RDFs for the solvent oxygen around the hydronium oxygen are shown in Fig. 4 where the result given by the ab initio molecular dynamics (AIMD) simulation 24 at the same thermodynamics conditions (ρ = 1.0 g/cm 3 , T = 300 K) is also provided as a reference. We note that the AIMD simulation in the graph was also yielded with the same exchange-correlation functional (BLYP) as our DFT calculation in the QM region. In Fig. 4  unit of the elementary charge. Thus, the QM water molecules in the first solvation shell are more polarized than the SPC/E water, though it is not straightforward to compare the QM charges with empirical MM partial charges. Therefore, it is possible that the water molecules relevant to the formation of the second peak will be attracted more strongly to the first solvation shell, which leads the decrease in the O-O distance and the enhancement of the distribution at the second peak. Of course, the water molecules in the second solvation shell are not treated quantum mechanically in our QM/MM-BCC simulation in contrast to the AIMD simulation. However, the SPC/E model of water will work adequately outside the first solvation shell since the effect of the ionic solute might be rather weakened. It should also be noted that the first peak in the RDF of QM/MM-BCC is narrower and its height is more emphasized than those of AIMD. We speculate that the distributions will be more to the bias potential is identical to the MM solvent. In the BCC method various types of constraint potential will be used and various forms of the shell Ω can be adopted according to the shape of the solute molecule of interest.
As benchmark tests we applied the BCC method to two kinds of water solutions. First we considered a solute water molecule embedded in an aqueous solution represented with a fully classical force field. We, then, computed the O-O RDFs by means of the BCC procedures for two different sizes of Ω (r c = 3.5 and 7.5Å). It was demonstrated that the RDFs given by the BCC method were in good agreements with those for the references. It should be kept in mind, however, that the choice of the number N A of the solvent molecules for a given size of the shell will somewhat affect the BCC results. Our recommendation is to take the appropriately smaller value than the average number N av of the solvent molecules inside the shell. More explicitly, N A should be determined so that it allows the deviation of the instantaneous number N in in the shell from N av . The number N av and the deviation can be readily obtained through a preliminary simulation. Further, we combined the QM/MM simulation with the BCC method (QM/MM-BCC) and applied it to the hydration of a hydronium ion. The O-O RDF for the solvent around the solute was computed using our method, which shows a fairy good agreement with a reference RDF given by an AIMD simulation.
Thus, it was shown that the BCC approach is simple yet effective and robust in computing a statistical property of a QM/MM system subjected to a constrained force. Our next issue as an extension of QM/MM-BCC is to compute the solvation free energy of an anionic QM solute in an aqueous solution. The formulation of the method and the test calculations are now proceeding. The results will be reported in forthcoming issues.