Intermediate electrostatic field for the elongation method

A simple way to improve the accuracy of the fragmentation methods is proposed. The formalism was applied to the elongation (ELG) method at restricted open-shell Hartree-Fock (ROHF) level of theory. The α-helix conformer of polyglycine was taken as a model system. The modified ELG method includes a simplified electrostatic field resulting from point-charge distribution of the system’s environment. In this way the long-distance polarization is approximately taken into account. The field attenuates during the ELG process to eventually disappear when the final structure is reached. The point-charge distributions for each ELG step are obtained from charge sensitivity analysis (CSA) in force-field atoms resolution. The presence of the intermediate field improves the accuracy of ELG calculations. The errors in total energy and its kinetic and potential contributions are reduced by at least one-order of magnitude. In addition the SCF convergence of ROHF scheme is improved.

In this paper we describe a simple way to improve the accuracy of the fragmentation methods. The self consistent field (SCF) calculations for selected molecular fragments are performed in the electrostatic field created by the point-charge distribution of the fragment's complementary part. The formalism is applied to elongation method where it is an intermediate construction. The field diminishes in each step of the elongation procedure and finally disappears after the electronic structure of the whole system is obtained. This intermediate electrostatic field introduces the long-distance polarization into ELG calculations and is a step beyond the nearsightedness approximation. There have been some attempts to include the polarization effects in charged polymers, however, they limited the fragment's complementary part to charged groups only. Here, we propose a uniform treatment of the fragment's complementary part. Point charge distribution was also used by Exner and Mezey [45] to diminish the distance criterion in the adjustable density matrix assembler method.
The paper is organized as follows. First, we briefly describe the ELG methods and charge sensitivity analysis (CSA) used to derive the point-charge distributions. Next, we give information connected with computational details. Then, the performance of the method is shown. The analysis is focused on the accuracy. Finally, the conclusions are given along with the future prospects.

Elongation and elongation cutoff technique
Elongation method mimics the polymerization/ copolymerization reaction mechanism [35,36]. The electronic structure of the whole molecular system (M) is synthesized by enlarging a so called starting cluster (M 1 ). The SCF calculations performed on M 1 initialize the ELG process. This step is complemented with molecular orbital (MO) localization procedure. The canonical M two fragments: A 1 and B 1 [37]. Next, the size of the system is enlarged by adding a new molecular fragment C 1 . The localized molecular orbitals (LMOs) assigned to fragment A 1 , which is far away from the chain propagation center, are frozen. The LMOs of B 1 and guess MOs of C 1 constitute the variation space S 1 . Then, SCF calculations on S 1 are performed. Such propagation procedure is continued, step by step, until M is built. Every SCF step of chain propagation is followed by molecular orbital localization. Canonical MOs are localized onto A i and B i fragments. The LMOs assigned to A i are excluded from variation space (they are frozen in the ELG process). The LMOs assigned to B i together with MOs of C i constitute the active space (S i ). The whole ELG procedure can be summarized as follows: Equations 1a, 1b, and 1c correspond to initialization, chain propagation, and chain termination, respectively. Notice that there is no need to localize MO in the final step (Eq. 1c).
In the ELG scheme the size of variation space is almost constant. One can eventually take advantage of the sparseness in the LMO representation [39,40,52,53]. In the limit of perfect localization, LMOs assigned to the frozen fragment (C A i ) have no tails in the active fragment and vice versa, LMOs assigned to the active fragment (C S i ) have no tails in the frozen one. Therefore, instead of constructing the Fock or Kohn-Sham matrix of M i-1 (F AO M i−1 M i−1 ), the F AO S i S i block is built. Such elongation cutoff scheme (ELG/C) operates within a low-dimensional subspace, and avoids the known bottleneck of the SCF calculations, i.e., diagonalization. The lower indices in Eq. (2) are introduced to indicate dimensions and the block 0 A i S i is filled with zeros. In addition, the number of two-electron repulsion integrals is substantially reduced as long as the total energy of M i is not needed. It was demonstrated that ELG/C scheme is linear in CPU time at HF [52] and KS [53] levels of theory for a linear or quasi-linear polymer. The ELG scheme can be generalized to three-dimensional (3D) systems. In such a generalized ELG scheme the frozen LMOs of a given fragment can reenter the variation space [54][55][56].
In the ELG (ELG/C) scheme the molecular fragments are not treated equivalently. Namely, the starting cluster does not know its "future" while S n-1 possesses the whole knowledge of M. Such way of building M is consistent with the nearsightedness approximation [28] which is a common assumption in the fragmentation methods. However, it is also a source of errors. To reduce this error, the ELG (ELG/C) process can be performed in an approximate field arising from complementary part of M i , for example its point charge distribution. The ith subsystem's complementary part is denoted by M i . By enlarging the system size (i→n) the complementary part becomes smaller and finally (i=n) disappears (M n is a zero set). This electrostatic field is an intermediate construction, therefore, the final energy of M should be greater than the reference HF (KS) one.

Charge sensitivity analysis in force-field atoms resolution
The intermediate point-charge distribution can be easily obtained from charge sensitivity analysis [57]. CSA is based on second-order Taylor expansion of the system's energy E M with respect to atomic charges. The CSA formalism in global resolution (without constraints on charge flow) can be summarized in a single matrix equation [58]: where η={η ij =∂ 2 E M /∂q i ∂q j } is the hardness matrix, q is the total charge and χ=∂E M /∂q is the global electronegativity [59, 60] of a molecular system M composed of N atoms. Vectors q=(q 1 ,q 2 ,…q N ) T and χ=(χ 1 * ,χ 2 * ,…χ N * ) T group the atomic charges and electronegativities, respectively. The first equation in (3) is a closure relation: The remaining equations are the electronegativity equalization equations (χ 1 =χ 2 = …=χ N =χ) [61]. The charge distribution inside M can be obtained by inverting Eq. (3): Here, β={β ij =∂ 2 E M /∂v i ∂v j =−∂q j /∂v i =−∂q i /∂v j } is the polarization matrix (linear response matrix). The vector v=(v 1 , v 2 ,…v N ) denotes the external potential due to nuclei. The remaining undefined quantities are the global hardness [62,63] η=∂ 2 E M /∂q 2 and the Fukui function (FF) [63] vector f={f i =(∂q i /∂q) v =−(∂μ/∂v i ) q }. The hardness matrix η and the vector of atomic electronegativities χ are the base parameters of CSA. In the force-field atoms resolution they depend on the atomic number, hybridization and local chemical environment of atoms constituting the system. We have recently derived these parameters for different population analyses [64,65].

Computational details
All ELG (ELG/C) calculations were performed at restricted open-shell Hartree-Fock (ROHF) level of theory using GAMESS package [66,67]. Three different basis sets, namely: STO-3G, 6-31G, and 6-31G(d), were applied. The alphahelix conformer of polyglycine was taken as a model system. The starting cluster was built of 15 amino acid units. In each step of the elongation eight units were frozen and another eight units were added to the system. This ELG propagation scheme is denoted as 15/8. The elongation process was terminated for 55 glycine units. In the case of ELG/C calculations, the cutoff procedure was initialized for a system made up of 31 units (15/8:31). The geometry of alpha-helix was the same as in our earlier paper [52]. We have chosen this conformer since, for a given partitioning scheme, the error in its total energy was bigger than for other conformers (C5, C7, and 3 10helix).
The point charge distribution for different population analyses was obtained from CSA calculations. Five different charge distributions were used, namely, Bader (B) [68], Hirshfeld (H) [69], Mülliken (M) [70], Natural (N) [71], and Voronoi (V) [72] population analyses. The ELG calculations were carried out for each charge distribution. The errors in total (tot), kinetic (kin), and potential (pot) energies were defined with respect to the conventional (supermolecule) ROHF energies: x ¼ tot; kin; pot ð Þ ð 7Þ By definition ΔE M =ΔE tot is greater than zero since ELG (ELG/C) is a variational method. There are no such restrictions on its components, therefore, ΔE kin (ΔE pot ) may be either negative or positive. We have chosen ROHF calculations since the propagation at RHF level of theory would require saturation of "broken" bonds. This means that each intermediate subsystem (15,23,31,39, and 47 units) would have to be saturated by a hydrogen atom. To avoid perturbation created by this artificial hydrogen atom ROHF scheme was selected. One can expect that the intermediate electrostatic field should stabilize the system and should improve the SCF convergence. The proposed modification of ELG (ELG/ C) scheme has no influence on CPU time since one-electron integrals are computed only once, before the SCF process.

Results and discussion
The errors introduced by ELG and modified ELG schemes are plotted in Fig. 1. Parts a, b, and c correspond to STO-3G, 6-31G, and 6-31G(d) basis sets, respectively. The error in the total energy is always positive. The ELG and modified ELG methods are variational, therefore the calculated energies are higher than the reference ROHF energies. All energies are computed for the final system built up of 55 units. The intermediate energies cannot be directly compared with the reference ROHF energies since they have quite a different meaning. The first error bar in each figure corresponds to standard ELG calculations with q M i ¼ 0 (no field). The remaining error bars correspond to different population analyses employed in the modified ELG procedure. The error in ELG calculation depends on the basis set and increases with its size. For 6-31G(d) basis set it is about two-orders of magnitude greater than for the minimal basis set. The presence of the intermediate field arising from charges from every population analysis reduces the error in the total energy. The reproduction of conventional ROHF energies is improved by about one order of magnitude for STO-3G (Fig. 1a) and 6-31G (Fig. 1b) basis sets. The improvement for 6-31G(d) basis is the most pronounced (Fig. 1c). The error is reduced by three orders of magnitude for B, M, and N population analyses. For all basis sets the modified ELG scheme with V and H charges works slightly worse than with B, M, and N charges. The changeability in total energy with the basis set size is the smallest for B and N charges. The errors in total energies are equal to 0.01, 0.02, and 0.02 kcal mol -1 for STO, 6-31G, and 6-31G(d) basis sets, respectively. Therefore, by taking the long-distance polarization into account, even in the simplified way, the reproduction of the reference ROHF energy is much better and the error does not exceed 1 kcal mol -1 . It should also be mentioned that, although Fig. 1 corresponds to ELG and modified ELG schemes, they present the errors in the ELG/C and modified ELG/C schemes. The reason for that is that the cutoff error (E tot ELG − E tot ELG/C ) is at least two orders of magnitude smaller than the elongation error (E tot ELG −E tot ROHF ). To understand the improvement in reproduction of the system's total energy we have decomposed E tot ELG into its potential (E pot ELG ) and kinetic (E kin ELG ) components. The errors with respect to conventional potential (E pot ROHF ) and kinetic (E kin ROHF ) energies are illustrated in Fig. 2. The black error bars correspond to ΔE kin while white error bars to ΔE pot . Again, the first error bar in each figure corresponds to standard ELG scheme. As it was mentioned in Computational details, the errors in potential and kinetic energies may be either positive or negative. It is clearly seen in Fig. 2a. Depending on the population analysis, the error in kinetic energy is positive (B and N) or negative (H, M and V). Regardless of the basis set and the population analysis, the error in potential energy always has the opposite sign to ΔE kin . Its magnitude is almost the same as that of ΔE pot . This error cancelation causes the accuracy in the total energy to be one order of magnitude greater than in its potential and kinetic components. Based on the virial theorem, a different ratio ΔE pot /ΔE kin should be expected. However, one should remember that the virial theorem is exact for the true ground-state wave function. The approximate wave function fulfils it only approximately. For more extended basis sets [6-31G and 6-31G(d)], ΔE pot for modified ELG scheme is negative and ΔE kin is positive. In the case of standard ELG scheme, ΔE pot is negative for 6-31G and positive for 6-31G(d) basis sets, respectively. It is clear from the figures that the intermediate electrostatic field improves the reproduction of kinetic and potential energies. Such improvement is evident for 6-31G(d) basis set. One should expect such behavior since slightly populated polarization functions on heavy atoms should be sensitive to the intermediate field used in the modified ELG calculations. All the results of modified ELG calculations at 6-31G(d) basis set show that the kinetic energy is overestimated due to the influence of point charges.
Finally, let us analyze the SCF convergence during ELG calculations. All intermediate subsystems are radicals (dublets). Only the whole system is closed-shell. The SCF convergence for radicals is worse than for closed-shell singlet states. The SCF convergence is illustrated in Fig. 3. Parts a, b, and c correspond to STO-3G, 6-31G, and 6-31G(d) bases, respectively. One can observe that, except for B, M, and N charges at STO-3G basis set calculations, the field stabilizes the radicals. It means that fewer SCF cycles are required to reach the same assumed accuracy in comparison to standard ELG scheme. The differences between ELG and modified ELG schemes are more pronounced for the starting cluster. For larger basis set, the M and N population analyses give the fastest convergence. In general, the natural population analysis (N) gives the best performance of calculations both in accuracy and convergence. At the end of the elongation process the number of iterations is the same for all curves. The system is then closed-shell and the number of iterations in the SCF cycle is small.

Conclusions
A simple modification of the elongation (elongation cutoff) method was proposed. Namely, the intermediate electrostatic field was introduced. The field is exerted by distributed monopoles located in the positions of atoms of the system's complementary part and it disappears in the final stage of the elongation calculations. Therefore, it does not violate the variational character of the ELG method. The modified ELG  scheme was tested for alpha-helix of polyglycine chain (55 glycine units). Several population analyses were applied. Charges were computed using CSA scheme, independently of the quantum chemical calculations. It was shown that the long-distance polarization, introduced by the field, improved the performance of the ELG method. The errors in the total, kinetic, and potential energies were reduced by at least one order of magnitude. The natural, Bader, and Mülliken population analyses gave the best agreement with the reference, conventional ROHF energies for the largest basis set. The proposed method can be easily adopted to other fragmentation techniques.
The modified ELG method improved the convergence of the SCF process at ROHF level of theory. The temporary electrostatic field stabilized the intermediate radical, therefore fewer cycles were involved during the SCF step. We plan to adopt the formalism at RHF and UHF levels of theory.