Efficient Geometry Minimization and Transition Structure Optimization

Subscriber access provided by Gothenburg University Library

Article

Efficient Geometry Minimization and Transition Structure Optimization Using Interpolated Potential Energy Surfaces and Iteratively Updated Hessians Jingjing Zheng, and Michael J. Frisch J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 18 Oct 2017 Downloaded from http://pubs.acs.org on October 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

revised to J. Chem. Theory Comput. October 16, 2017 Efficient Geometry Minimization and Transition Structure Optimization Using Interpolated Potential Energy Surfaces and Iteratively Updated Hessians Jingjing Zheng* and Michael J. Frisch Gaussian, Inc., Wallingford, Connecticut 06492 Abstract An efficient geometry optimization algorithm based on interpolated potential energy surfaces with iteratively updated Hessians is presented in this work. At each step of geometry optimization (including both minimization and transition structure search), an interpolated potential energy surface is properly constructed by using the previously calculated information (energies, gradients, and Hessians/updated Hessians), and Hessians of the two latest geometries are updated in an iterative manner. The optimized minimum or transition structure on the interpolated surface is used for the starting geometry of the next geometry optimization step. The cost of searching minimum or transition structure on the interpolated surface and iteratively updating Hessians is usually negligible compared with most electronic structure single gradient calculations. These interpolated potential energy surfaces are often better representations of the true potential energy surface in a broader range than a local quadratic approximation that is usually used in most geometry optimization algorithms. Tests on a series of large and floppy molecules and transition structures both in gas phase and in solutions shows that the new algorithm can significantly improve the optimization efficiency by using the iteratively updated Hessians and optimizations on interpolated surfaces. .

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

I.

Page 2 of 29

Introduction Geometry optimization takes an essential role in most electronic structure calculations.

Various geometry optimization algorithms have been implemented in many computational chemistry software packages. Geometry optimization in the context of electronic structure calculations usually refers to the search of the minimum-energy structure and transition structure. Several reviews about geometry optimizations are given in refs. 1,2,3,4,5,6,7,8. Improving efficiency is a main goal of developing geometry optimization methods. How to maximizing the usage of the already computed information on the potential energy surface is a key to improve the efficiency of optimization. In most algorithms, the step to the desired stationary point is determined by the gradient and the Hessian of the current point on the potential energy surface (PES). For example, the Newton or Newton-Raphson (NR) step gives the displacements ∆ x to the minimum under the quadratic approximation by ∆x = − H −1g

(1)

Where H is Hessian matrix and g is gradient vector computed at the current step. When Hessian is obtained in an approximate way, the Newton method is usually called quasi-Newton method. The NR method has enhanced stability by employing the techniques such as rational function optimization (RFO)9,10 and the trust radius model.1,4,11, 12, 13 One often used technique to reduce the number of optimization steps is the direct inversion in the iterative subspace (DIIS).14, 15 Geometry optimization by DIIS (GDIIS)16,17 or energy-represented DIIS (GEDIIS)18 method adopts a set of vector { R i } including information of geometry and gradient previously calculated and extrapolates/interpolates { R i } by minimizing the errors in the least-square manner. In a simultaneous optimization method19 for 2 ACS Paragon Plus Environment

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


both molecular geometry and electronic wave function, a constant Hessian is included in the vector { R i }. Usually GDIIS or GEDIIS for minimization converges faster than quasi-Newton methods at the region near minimum, e.g. where root mean square (RMS) of gradient is smaller than 2.5 ×10−3. Therefore a hybrid of quasi-Newton and GDIIS (GEDIIS) methods can be more efficient.18 The RFO and GDIIS (GEDIIS) methods use the previously obtained PES information in some ways. For example, Hessians are usually updated by using the previous gradients in these methods; the DIIS and GEDIIS methods use some previously calculated geometry, energy, and gradient to do linear interpolation/extrapolation. In this paper the new idea about recycling the already calculated points is to use those points to fit an analytic potential energy surface to represent the true potential energy surface and to guide the geometry optimizations. And the updated Hessians are improved in an iterative and balanced manner. The purpose of fitting potential energy is different with that used in the dynamics where the fitted surface should be as precise as possible for the dynamical interested region. Usually the number of configurations is very limited during the geometry optimization compared to very large number of configurations used in the regular fitting potential energy surface process, and also these configurations do not cover the desired minimum or transition structure. Therefore our purpose here to fitting a potential energy surface is to obtain a better representation for a larger range of potential energy surface than that represented by a local harmonic approximation. The fitted potential energy surface is reconstructed at each geometry optimization step. In this work, we present an efficient method for both geometry minimization and transition state search using properly constructed interpolated potential energy surfaces with the improved updated Hessians.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

II.

Page 4 of 29

Method II.A. Constructing interpolated surfaces. At each geometry optimization step, we

construct an interpolated potential energy surface using a few previous consecutively calculated points (geometries, energies, gradients, and Hessians). Then we perform optimization on the interpolated potential energy surface to locate the desired stationary point (minimum or transition structure). The displacement taken on the true potential energy surface is the difference between the current geometry and the optimized geometry on the interpolated surface. Note that a complete geometry optimization process on the interpolated surface has negligible computational cost compared with an electronic structure calculation for an energy and gradient. Let’s call the optimization on an interpolated surface as micro-optimization in this paper to distinguish the optimization on the true potential energy surface. The interpolated potential energy surface V int is constructed as N

V int = ∑ WiVi

(2)

1 Vi = V0 (R i ) + giT (R − R i ) + (R − R i ) T H i (R − R i ) 2

(3)

i =1

where N is the number of previously calculated points used in the interpolation, Wi is the normalized weighting function for the point i, the potential V0 (R i ) is the potential energy at geometry R i , g i is gradient, and H i is Hessian. Note that the R i can be any coordinates and is not necessary to be the same as the coordinates used for geometry optimization; in this work,


Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


redundant internal coordinates are used. The normalized weighting function Wi is given by normalization of the un-normalized weighting function wi

Wi =

wi N

(4)

∑wj j =1

There are numerous ways to weight each point. In this work, we categorize the interpolations according to the weights into three types: energy-gradient weighted interpolation (EGWI), distance weighted interpolation (DWI), and combined distance-energy-gradient weighted interpolation (DEGWI). Energy-Gradient Weighted Interpolation (EGWI). Here we present three weighting functions for EGWI; they are exponential function, rational function, and hyperbolic tangent function, respectively. Exponential weighting function. The un-normalized exponential weighting function wi is

wiEG − Exp = exp[− ∆ui ]  G RMS (g )  V (R ) − V i  i min ∆u i = C ln  G RMS  (Vmax − Vmin )  min 

(5a)

(5b)

for geometry minimization, where C is a constant that controls the decay rate of the weight, RMS , and G RMS (g i ) is the root mean square of the components in gradient g i , Vmin , Vmax , G min

RMS are minimum and maximum values of the potential V (R i ) and G RMS (g i ) , respectively. G max

For transition structure search, the function ∆u i is only the function of gradients



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

 G RMS (g )  G RMS (g i ) − G RMS i  min ∆ui = C ln  G RMS  G RMS − G RMS max min min  

Page 6 of 29

(6)

The constant C for exponential weighting function is parameterized to be 5.5 using the test suites in Table 1 and 3. Rational weighting function. The un-normalized rational weighting function using rational function is

wiEG − RF =

1 1 + ∆ui

(7)

The constant C in Eqs. 5b and 6 for rational weighting function is parameterized to be 40.0 using the test suite in Table 1 and 3. Hyperbolic tangent weighting function. The un-normalized rational weighting function using hyperbolic tangent (hypertangent) function is wiEG − HT = 1 − tanh[ ∆ui ]

(8)

The constant C in Eqs. 5b and 6 is parameterized to be 2.5 for hyperbolic tangent weighting function using the test suites in Table 1 and 3. Distance Weighted Interpolation (DWI). The distanced-weighted interpolation is also known as Shepard interpolation20,21 that has been successfully applied to fit potential energy surfaces that are often used for molecular dynamics. 22,23,24,25,26,27,28,29 The weighting function used in Shepard interpolation is a function of geometry, and therefore the derivatives of the interpolated potential also become more complicated.


Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The un-normalized weighting function used for minimization and transition structure search have the form

wiD =

1 x − xi

2p

(9)

Where p is an integer and is chosen as 2 in our calculations. Distance-Energy-Gradient Weighted Interpolation (DEGWI). Since we already have the interpolation schemes based on either distance or energy-gradient associated functions, an alternative way to weight the interpolation points is to use all the information of distance, energy, and distance. The DEGWI un-normalized weighting functions are taken as the product of DWI weighting function and one of the EGWI weighting function, such as

wiDEG = wiD wiEG − Exp

(10)

wiDEG = wiD wiEG − RF

(11)

wiDEG = wiD wiEG − HT

(12)

The local quadratic approximation of potential energy surface is only often valid for a very small region. These interpolated surfaces expand the local quadratic approximation to a broader region by including gradients and Hessians of surrounding points. II.B. Iteratively Updated Hessians. In geometry optimization, it is well known that the accurate Hessian is critical to the efficiency of geometry optimization. However explicitly calculating Hessian in each step is usually too expensive therefore Hessians are usually obtained 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

by using various updating schemes. In this study, Hessians are updated with first derivatives using the Broyden-Fletcher-Goldfarb-Shanno (BFGS)30, 31, 32, 33 scheme for minimizations and using the Bofill’s weighted scheme34 for transition structure optimizations. Here updating Hessian means that a Hessian of a target geometry is calculated using the analytic gradients of that given geometry and one or more reference geometries. When updating Hessian using previous points and the target point, the updated Hessian of the target geometry is calculated as H tar = H ref + ∆H

(13)

The ∆H is obtained by a Hessian update scheme, e.g. BFGS or Bofill’s scheme. All the Hessian update schemes fulfill the Newton condition, ∆g = H tar ∆x

(14)

where ∆g = g(x tar ) − g (x ref ) and ∆ x = x tar − x ref . Note that the ∆H is a function of coordinates, gradients, and Hessians of a reference point in all the updating schemes. We also note that when a new geometry is calculated during the optimization, the new geometry can be a reference point for updating the Hessians of previously calculated points, and these improved Hessian of the last point including the new geometry information then become the reference Hessian of the new geometry. For computational efficiency, we only go one step back to update Hessian using the new geometry as a reference. When N reference points are available for a target point, we apply the eq 13 in a sequential manner according to the distances between the reference point and target point. Here we point out that the final updated H tar depends on the sequence of applying all the reference points when using the BFGS and Bofill’s methods. If the H iref is the updated Hessian after applying the (i-1)-th point, according to .eq. 13 it is


Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


H iref = H iref −1 + ∆H i −1

(15)

Note that ∆H i −1 is also a function of H iref −1 . When the sequence of the points is changed, each H iref would be changed, and eventually the updated H tar would be different. This conclusion is

numerically proofed by the fact that the number of optimization steps is changed when the sequence of the points is reversed in the Hessian updates (see Section IV for details). For the same reason mentioned above, when the same set of points are applied to update Hessian multiple times, the final updated H tar is different when they are applied only once. Therefore, we repeat the whole sequential process N times (and we set 50 times at maximum in our implementation). That is to say, the eq. 13 is applied N2 times for a given point. We call this Hessian updating scheme as iterative updating (IU) in the rest of paper. When applying the GIS optimization, we always use the IU for updating Hessians, so that the points used in the interpolations have the same quality Hessians. It is crucial because using a lower quality Hessian in the interpolation could downgrade interpolated surface than the local harmonic approximation. Applying the iterative updating scheme will increase the computational cost, but this cost of updating Hessian is still much less than a single point gradient calculation using an electronic structure method. The micro-optimization on the interpolated surface is carried out either by the RFO method or by a combination linear search and the RFO method.35 Since an analytic potential energy surface is used, the convergence criteria for micro-optimization is set to very tight condition, i.e. RMS internal force is smaller than 1 ×10−8 au and the RMS displacement is smaller than 1×10−5 au. We call this geometry optimization on interpolated surface GIS



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

The geometries of the first few steps during optimization are usually far away from the desired stationary point, so it would not be helpful to construct an interpolated surface to guide the optimization using the information from the geometries too far away from the stationary point. Therefore we adopted a hybrid scheme to perform minimization and transition structure optimization. This hybrid scheme starts optimization using RFO method and switches to the GIS method when RMS force in internal coordinates is less than 5×10−3 au. III. Computations and Results All calculations are performed using the development version of the Gaussian series of programs.36 The convergence criteria is met when the RMS force is less than 1×10−5 au, the RMS geometry displacement is less than 4×10−5 au, the maximum component of the force vector is less than 1.5×10−5 au, and the maximum component of geometry displacement is less than 6×10−5 au; we call this convergence criteria as tight convergence. Note that the above force and displacement are in Cartesian coordinates although internal coordinates are used in in optimization iterations. The default convergence criteria in Gaussian program is also used to access the convergence behavior, in which the optimization is considered as converged when the RMS force is less than 3×10−4 au, the RMS geometry displacement is less than 1.2×10−3 au, the maximum component of the force vector is less than 4.5×10−4 au, and the maximum component of the geometry displacement is less than 1.8×10−3 au. The redundant internal coordinates are used in the optimizations.37,38,39,40,41,42 For the purpose of comparison, we have also carried out geometry optimizations for the same set of molecules and transition structures using the GEDIIS method for minimization and the RFO method for transition structure optimization, which are the default algorithms for 10 ACS Paragon Plus Environment

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


geometry optimization implemented in the Gaussian 16 program. The mixed scheme is the same as that in the ref. 18, that is (1) the geometry optimization begins with the RFO method to use the Hessian-based quadrature to avoid shallow potential wells; (2) the optimization scheme switches to the GDIIS when RMS force in internal coordinate is smaller than 2.5×10−3 au; (3) when optimization is switched to GDIIS scheme and the current point energy is not the lowest energy point among the points used in the GDIIS, the GDIIS is replaced by the GEDIIS method. We simply call this mixed method as GEDIIS method through the rest of the paper.

III.A. Minimization. A set of large molecules are used to test the performance of geometry minimization either on AM143,44 potential energy surfaces and/or PM645 potential energy surfaces. These molecules are, α-tocopherol, γ-cyclodextrins, taxol, a series of alanine peptides, i.e. (Ala)5, (Ala)9, (Ala)10, and (Ala)25, and various sizes of water clusters, i.e. (H2O)10, (H2O)16, and (H2O)25. The sizes of these systems vary from 30 atoms to 259 atoms. The numbers of steps required to converge in each optimization are shown in Tables 1. In Tables 1 and 4, all GIS calculations employ two points (the current step and the point from one step before the current one) to construct an interpolated surface in each geometry optimization step. We also performed optimization using GEDIIS method with the iterative updating Hessians in order to see the effect of the IU scheme (GEDIIS+IU). Note that in the GEDIIS+IU optimization, only the Hessian of the current point is used for the quadratic step. The linear extrapolation (LE) in Berny’s algorithm35 also uses two points for extrapolation, but only use the energies and gradients. Therefore we also perform the geometry minimizations using combination of linear extrapolation and RFO technique for the same of molecules, which is denoted as RFO+LE. The results are shown in Table 1. 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

To illustrating the optimization processes, we plot the energies, RMS force in internal coordinates, and RMS of displacements along the optimization path for taxol (on AM1 surface) cases in Fig. 1. The calculations in these plots employ the exponential weighting function and use two points in the interpolations.

III.B. Transition Structure. Six transition structures from Ref. 46 and seven transition structures from Ref. 47 are used to compare the GIS and RFO optimization methods. Table 2 summarizes the reactions and the methods for generating initial structure, initial Hessians, and solvent and solvation models. The APFD48 density functional method with 6-31+G(d,p) basis set are used for the transition structures in Ref 46 and PM6 method is used for transition structures from Ref.47. The numbers of optimization steps for both tight and default convergence criteria are listed in the Table 3. The RFO method with IU scheme was also used to perform the optimization for comparison purpose. The energies, RMS internal forces, and RMS of displacements along the optimization paths are plotted in Fig. 1 for the Grignard case.

IV. Discussion In order to understand how the GIS method works in a simple picture, a one-dimension (1D) anharmonic surface is used to test this idea as well as the Newton-Raphson method. The 1D surface is

y = x + x3 + x4

(16)

This surface is shown in Figure 2 and it has a quite flat region of [-1.5, 0.5] that is typically difficult for geometry optimization. And also the large anharmonicity of this surface could make


Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the movement by the NR method based on quadratic approximation too large or too small. We set the initial value of x to be 1.9, and applied both NR method and NR-GIS hybrid method to search the minimum of this curve. In applying NR method, the gradient and second derivative of the function y is calculated explicitly at each step. In applying GIS method, when the actual analytic gradient is smaller than 1.4 the GIS method is used otherwise the NR method is used. Two previous points are used for constructing an interpolated surface. The un-normalized weighting function is simplified as

[

wi = exp − 2( y − y min ) 2

]

(17)

where i = 1 or 2, ymin is the minimum value between y1 and y 2 . In applying the NR method, the displacement ∆x = − y ′′ / y ′ , while in the NR-GIS method, the derivatives of y are calculated on the interpolated surface. Figure 3 shows the displacement ∆x , value of y, first and second derivative of y along the optimization path. The NR-GIS method convergence faster than NR method, and the values of y and y ' vary much smoother along optimization path than the NR method. At the step 4, the function has a very small negative value (−0.65) of second derivative so that it produce a quite large step for the wrong direction, which cause a big bump at step 5. However, the interpolated surface at step 4 has the positive second derivative (0.70) by taking the weighted average of y ' ' at step 4 and 3, therefore the curvature of the interpolated surface has better representation for the whole curve and make the movement go to the right direction. Figure 3 also shows the interpolated curve (green dash line) at step 4 and the curve calculated by the quadratic approximation (red dotted dash line) at step 4. These two curves shows that the interpolated



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 29

surface can be considered more global than the local quadratic approximation by taking the information from the previous steps. This 1D model shows how an interpolated surface in optimization woks in a very straightforward manner. It also shows how the efficiency of optimization could be largely improved by using properly constructed interpolated surface. In general, Table 1 shows that the GIS method outperforms the GEDIIS method and RFO+IE method for these complex molecules under both tight and default convergence criteria. When the IU scheme is combined the GEDIIS method, it is also quite obvious that the iterative updating scheme does improve the Hessians so that the efficiency of optimization is improved compared with the GEDIIS method. For quite a few cases, the GIS can significantly improves the optimization efficiency, e.g. (Ala)9 on AM1 surface and (Ala)25. Especially at the near convergence regions, one needs more precise information of PES in order to reach the tight convergence. Table 1 shows that the performance of GIS method leading to the tight convergence is more stable than that GEDIIS and RFO+IE methods, which is due to the interpolated surface including the energies, gradients, and Hessians of the latest two points in a proper way. Figure 1 shows the example that GIS optimization leads energies, RMS forces, and RMS displacements going down more smoothly and optimization runs in a more controlled way. We also notice that different methods could lead optimization to different minima for a complex molecule. But neither of the method is designed to lead optimization to lower-energy minima. In Section II.B, we mentioned that the updated Hessian depends on the sequence of reference points. In our IU scheme, we apply the reference point according to their distance to the target point, i.e. from the farthest point to the closest point. To proof the dependence, we


Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reverse the sequence and find that the number of optimization steps to convergence usually changed. For example, (Ala)5 on AM1 surface needs five more steps to converge. The GIS method performs also better than the RFO method for almost all the cases of transition structure optimizations as shown in Table 3. For many cases, the improvement by the GIS method is dramatic and much of the improvement are benefited from the IU scheme. Figure 1 shows the GIS method gives a much smoother and stable optimization paths for the Grignard case. Table 4 shows that various weighting schemes have similar performance for most cases. So we conclude that the improved efficiency of the GIS method is not very sensitive to the weighting as long as the weights are in a reasonable range. All the tested weighting schemes seem reasonable, however, it is difficult to conclude one scheme is better than the other from these limited testing cases. So far, we have seen the improvements of the GIS method using two points in the interpolations. We also test the GIS optimizations for the same set of molecules and transition structures using three or more points in the interpolations. However, we did not observe systematic improvement by using more points. That is to say, sometimes more points in the interpolation can enhance the efficiency, but more point could also downgrade the efficiency for some cases if some noisy points are included. So how to choose points in the interpolation in a more intelligent way or how to construct a more robust weighting scheme still needs more thoughts and better understanding of the GIS method. In a regular potential energy surface fitting process, often a very large number of geometries, energies, and/or gradients (e.g. thousands or even much more) are used to fit a



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 29

surface covering all input geometries as close as possible to the true surface. However in geometry optimizations the number of available geometries and energies are often quite small and most geometries are far away from stationary point, therefore it is difficult to follow the regular surface fitting protocol in order to obtain an accurate surface that cover the minimum or transition structure region. Note that the purpose of improving the optimization efficiency is to find the stationary point in the least steps. In the GIS method, we try to obtain an analytic surface that represent better curvature of the local area so that this surface can lead the stepping to the right direction and then improves the efficiency of geometry optimization. This interpolated surface may not be close to the true surface when the interpolants are far away from stationary points, but it is reconstructed at each step and eventually converges to the true surface. V.

Conclusion In this article, we present a new geometry optimization method based on properly

constructed interpolated potential energy surfaces and using iterative updated Hessians. These interpolated potential energy surfaces are often better representation of the true potential energy surface in a broader range than a local quadratic approximation that is usually used in most geometry optimization algorithms. By having a better approximation of PES at each step, the GIS method improves the efficiency of geometry optimization for both minimization and transition structure search. In this work, we present three class of weighting scheme based on distance, energy, and/or gradient. The energy-gradient based weighting functions adopts the function form as exponential function, rational function, and hyperbolic tangent function. The results presented in this work show that the optimization efficiency is not very sensitive to the form of weighting functions and the parameters used in the weighting functions. There still could


Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


be more room to improve the efficiency of geometry optimization by constructing interpolated potential energy surfaces in more intelligent ways, which is a subject of future research.

Author Information Corresponding author *Email: [email protected]



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 18 of 29

Table 1 Numbers of geometry minimization step and energies for a set of molecules.

____________________________________________________________________________________________________________ Molecule

Default Convergence GEDIIS

GEDIIS+IUa

RFO+LEb

Tight Convergence GISc

GEDIIS

GEDIIS+IUa

Final Energy (hartree)

RFO+LEb

GISc

GEDIIS

GEDIIS+IUa

RFO+LEb

GISc

_________________________________________________________________________________________________________________________________ (Ala)25 AM1

87

81

85

76

117

125

511

84

-1.7793317 -1.7793317 -1.7793317 -1.7793317

(Ala)5 PM6

27

28

28

24

33

34

34

30

-0.4669919 -0.4669919 -0.4669919 -0.4669919

(Ala)5 AM1

40

37

37

35

47

47

47

40

-0.3975799 -0.3975799 -0.3975799 -0.3975799

(Ala)9 PM6

37

40

33

33

44

45

44

40

-0.8245757 -0.8245757 -0.8245757 -0.8245757

(Ala)9 AM1

63

41

41

32

74

50

45

39

-0.6741379 -0.6741379 -0.6741379 -0.6739792

α-tocopherol AM1

44

35

34

35

57

53

48

47

-0.2915054 -0.2915054 -0.2915054 -0.2915054

Decamethylzincocene AM1

32

27

24

31

35

33

32

38

γ-Cyclodoxtrins AM1

97

90

74

74

110

98

85

81

-3.0180088 -3.0185662 -3.0185662 -3.0185662

(Ala)10 Conformation 1 AM1

90

77

88

77

119

97

109

86

-0.7333443 -0.7333443 -0.7333443 -0.7331108


89

69

68

66

103

85

93

76

-0.7337492 -0.7337492 -0.7337492 -0.7337492

Taxol AM1

64

54

52

54

81

96

196

61

-0.6668621 -0.6668621 -0.6668621 -0.6668621

Taxol PM6

57

54

56

52

72

83

NCd

66

-0.8047769 -0.8047769

(H2O)10 AM1

31

29

31

31

58

47

50

46

-0.9792547 -0.9792547 -0.9792547 -0.9792547

(H2O)16 AM1

41

38

39

40

61

52

54

52

-1.5949848 -1.5949848 -1.5949848 -1.5949848

(H2O)21 AM1

102

101

97

96

141

122

117

115

-2.0900857 -2.0900857 -2.0900857 -2.0900857

0.0252448

0.0252448

0.0252448

0.0252448

-0.8047769

____________________________________________________________________________________________________________________________ aGEDIIS bEach

optimization step uses RFO for quadratic step with iteratively updated Hessian and linear extrapolation (LE) for linear step.

cUsing dNot

method with iteratively updating Hessians

the exponential weighting function and 2-points interpolation

converged.


Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Table 2 List of transition structures _____________________________________________________________________________________________________________________________ Reaction

Initial guess

Initial Hessian

Solvent

Model Chemistry

______________________________________________________________________________________________________________________________ TS Diss

Isopropylazide → dimethylamine + N2

QST2

Approx.

No

APFD/6-31+G(d,p)

TSInsert

SiH4…Si+ → SiH3SiH+

QST3

Approx.

No

APFD/6-31+G(d,p)

TS1BAC2

Methyl acetate + OH− + H2O → intermediate

QST3

Analytical

watera

APFD/6-31+G(d,p)

TS2BAC2

Intermediate → acetate ion + methanol

QST3

Analytical

watera

APFD/6-31+G(d,p)

N2 TSSOH

methyl acetate + OH− + H2O → methanol + acetyl ion

QST3

Analytical

watera

APFD/6-31+G(d,p)

TSSClN2

Chloroethane + acetate ion → ethyl acetate + Cl−

QST2

Analytical

DMSOb

APFD/6-31+G(d,p)

C5HT

Hydrogen transfer from C5 to C1 in 1,3-penadine

QST2

Estimated

No

PM6

Cope

Cope rearrangement of 1,4-hexadiene

QST3

Estimated

No

PM6

DACP2

Diel-Alder addition for cyclopentadiene and ethylene

QST3

Estimated

No

PM6

Ene

Ene reaction of ethylene and propene

QST2

Estimated

No

PM6

Grignard

Addition of phenyl magnesium bromide

QST3

Estimated

No

PM6

Oxirante

ring opening of 2-methyl-3-phenyloxirane

QST3

Esimated

No

PM6

Sulfolene

addition of sulfur dioxide to butadiene

QST2

Estimated

No

PM6

_________________________________________________________________________________________________________________________ a PCM model is used. b SMD model is used



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 29

Table 3 Numbers of stepsa for transition structure search using RFO and GIS methods. ___________________________________________________________________________________________ Default Convergence GEDIIS

GEDIIS+IUa

Tight Convergence GISb

GEDIIS

GEDIIS+IUa

GISb

____________________________________________________________________________________________ TS Diss

28

21

22

39

27

28

TSInsert

22

23

23

31

27

27

TS1BAC2

NCc

28

34

NCc

38

47

AC2 TSB 2

NCc

37

37

NCc

42

40

N2 TSSOH

60

37

63

92

NCc

78

TSSClN2

31

19

19

50

33

36

C5HT

15

19

15

29

21

20

Cope

33

33

32

46

35

34

DACP2

29

19

21

34

24

23

Ene

35

28

27

43

32

30

Grignard

66

43

37

81

53

48

Oxirane

NCc

NCc

96

NC

NCc

111

25

22

23

34

27

26

Sulfolene

______________________________________________________________________________ aGEDIIS bUsing cNot

method with iteratively updating Hessians

the exponential weighting function and 2-points interpolation

converged by reaching the maximum number of steps.


Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4 Numbers of optimization steps using GIS method with various weighting functions in the interpolations. Two points are used in the interpolations and tight convergence criteria is applied. _____________________________________________________________________________________________ DWI

Exp

RF

Hypt

DExp

DRF

DHypt

_____________________________________________________________________________________________ (Ala)25 AM1

96

84

87

85

93

88

87

(Ala)5 PM6

30

30

29

30

30

30

30

(Ala)5 AM1

41

40

40

40

46

42

43

(Ala)9 PM6

43

40

40

40

48

42

48

(Ala)9 AM1

42

39

40

39

38

38

42

α-tocopherol AM1

47

47

45

44

47

50

46

Decamethylzincocene AM1

38

38

38

38

35

35

35

γ-Cyclodoxtrins AM1

91

81

82

83

74

76

73


96

86

90

86

104

86

104


77

76

81

75

78

77

75

Taxol AM1

62

61

63

64

61

84

63

Taxol PM6

63

66

63

67

64

66

69

(H2O)10 AM1

47

46

47

47

45

45

45

(H2O)16 AM1

52

52

52

52

52

52

52

(H2O)21 AM1

112

115

114

115

114

112

114

______________________________________________________________________________



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. Energies and RMS internal forces along optimization paths for taxol and transition structure

Page 22 of 29

TSSClN2 .

SN2 Taxol is optimized on the AM1 surface and transition structure TSCl is optimized on the PM6 surface. The exponential weighting function with two points is used in constructing interpolated surface in RFO-GIS method.


Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2 The black solid line is a 1D anharmonic function by eq 12; the green dash line is the interpolated curve at step 4 and the red dotted dash line is the quadratic curve at step 4.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29

Figure 3 The displacement ∆x , value of y, first derivative y ' , and second derivative y ' ' along the optimization path. Note that the shown y and y ' for NR-GIS method are the actual values instead of the interpolated values, while y ' ' is the second derivative of the interpolated surface whenever the GIS is applied in the hybrid method. However, we do use y ' and y ' ' on the interpolated surface to generate the displacement ∆x when the GIS method is applied.


Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References

1

Schlegel, H. B. Geometry Optimization. WIREs Comput. Mol. Sci. 2011, 1, 790-806.

2

Hratchian, H. P.; Schlegel H. B. Finding Minima, Transition States, and Following Reaction Pathways on Ab Initio Potential Energy Surfaces. In Theory and Applications of Computational Chemistry: the First Forty Years; Dykstra, C. E.; Frenking G.; Kim, K. S.; Scuseria, G. E., Ed.; Elsevier, New York, 2005, pp 195-249.

3

Wales, D. J. Energy Landscape; Cambridge University Press, 2003.

4

Schlegel, H. B. Geometry optimization on potential energy surfaces. In Modern Electronic Structure Theory; Yarkony D. R., Ed.; World Scientific Publishing, Singapore, 1995, pp 459-500.

5

Liotard, D. A. Algorithmic tools in the study of semiempirical potential surfaces. Int. J. Quantum Chem. 1992, 44, 723-741.

6

Schlegel, H. B. Some practical suggestions for optimizing geometries and locating transition states. In New Theoretical Concepts for Understanding Organic Reactions; Bertrán, J., Ed.; Kluwer Academic, The Netherlands, 1989, 33-53.

7

Head, J. D.; Weiner, B.; Zerner, M. C. A survey of optimization procedures for stable structures and transition states. Int. J. Quantum Chem. 1988, 33, 177-186.

8

Schlegel, H. B. Optimization of equilibrium geometries and transition structures. Avd. Chem. Phys. 1987, 67, 249-286.

9

Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. Search for stationary points on surfaces. J. Phys. Chem. 1985, 89, 52-57.

10 Simons, J.; Nichols, J. Strategies for walking on potential energy surfaces using local quadratic approximations. Int. J. Quantum Chem. 1990, 24, 263-276. 11 Fletcher, R. Practical Methods of Optimizations; Wiley, Chichester, UK, 1981. 12 Murray, W.; Wright, M. H. Practical Optimization; Academic, New York, 1981. 13 Nonlinear Optimization; Powell, M. J. D., Ed.; Academic, New York, 1982. 14 Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chem. Phys. Lett. 1980, 73, 393-398. 15 Pulay, P. Improved SCF Convergence Acceleration. J. Comput. Chem. 1982, 3, 556-560.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 29

16 Csaszar, P.; Pulay, P. Geometry optimization by direct inversion in the iterative subspace. J. Mol. Struct. 1984, 114, 31-34. 17 Farkas, Ö.; Schlegel, H. B. Methods for optimizing large molecules. Part III. An improved algorithm for geometry optimization using direct inversion in the iterative subspace (GDIIS). Phys. Chem. Chem. Phys. 2002, 4, 11-15. 18 Li, X.; Frisch, M. J. Energy-Represented Direct Inversion in the Iterative Subspace within a Hybrid Geometry Optimization Method. J. Chem. Theory Comput. 2006, 2, 835-839. 19 Moss, C. L.; Li, X. First order simultaneous optimization of molecular geometry and electronic wave function. J. Chem. Phys. 2008, 129, 114102/1-6. 20 Ischtwan, J.; Collins, M. A. Molecular potential energy surfaces by interpolation. J .Chem. Phys. 1994, 100, 8080-8088. 21 Collins, M. A. Molecular potential-energy surfaces for chemical reaction dynamics. Theor. Chem. Acc. 2002, 108, 313-24. 22 Bettens, R. P. A.; Collins, M. A. Potential energy surfaces and dynamics for the reactions between C(3P) and H 3+ ( 1 A1' ) J. Chem. Phys. 1998, 108, 2424-2433. 23 Bettens, R. P. A.; Collins, M. A. Interpolated potential energy surface and dynamics for the reactions between N(4S) and H 3+ ( 1 A1' ) J. Chem. Phys. 1998, 109, 9728-9736. 24 Collins, M. A.; Zhang, D. H. Application of interpolated potential energy surfaces to quantum reactive scattering. J. Chem. Phys. 1999, 111, 9924-9931. 25 Fuller, R. O.; Bettens, R. P. A.; Collins, M. A. Interpolated potential-energy surface and reaction dynamics for BH++H2. J. Chem. Phys. 2001, 114, 10711-10716. 26 Kim, Y.; Corchado, J. C.; Villà, J.; Xing, J.; Truhlar, D. G. Multiconfiguration molecular mechanics algorithm for potential energy surfaces of chemical reactions. J. Chem. Phys. 2000, 112, 2718-2735. 27 Song, K.; Collins, M. A. A classical trajectory study of sym-triazine photodissociation on an interpolated potential energy surface. Chem. Phys. Lett. 2001, 335, 481-488. 28 Hratchian, H. P.; Schlegel, H. B. Accurate reaction paths using a Hessian based predictor– corrector integrator. J. Chem. Phys. 2004, 120, 9918-9924.


Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


29 Hratchian, H. P.; Schlegel, H. B. Using Hessian Updating To Increase the Efficiency of a Hessian Based Predictor-Corrector Reaction Path Following Method. J. Chem. Theory Comput. 2005, 1, 61-69. 30 Broyden, C. G. The convergence of a class of double-rank minimization algorithms. I: General considerations. J. Inst. Math. Appl. 1970, 6, 76-90. 31 Fletcher, R. A new approach to variable metric algorithms. Comput. J. 1970, 13, 317-322. 32 Goldfarb, D. A family of variable-metric methods derived by variational means. Math. Comput. 1970, 24, 23-26. 33 Shanno, D. F. Conditioning of quasi-Newton methods for function minimization. Math. Comput. 1970, 24, 647-657. 34 Bofill, J. M. Updated Hessian matrix and the restricted step method for locating transition structures. J Comput. Chem. 1994, 15, 1-11. 35 Schlegel, H. B. Optimization of equilibrium geometries and transition structures. J. Comput. Chem. 1982, 3, 214-218. 36 Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B; Petrone, A.; Ortiz, J. V.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, F; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Development Version Rev. I10+, Gaussian, Inc., Wallingford CT, 2016. 37 Pulay, P.; Fogarasi, G.; Pang, F.; Boggs, J. E. Systematic ab initio gradient calculation of molecular geometries, force constants, and dipole moment derivatives. J. Am. Chem. Soc. 1979, 101, 2550-2560.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 29

38 Fogarasi, G.; Zhou, X.; Taylor, P.; Pulay, P. The calculation of ab initio molecular geometries: efficient optimization by natural internal coordinates and empirical correction by offset forces. J. Am. Chem. Soc. 1992, 114, 8191-8201. 39 Pulay, P.; Fogarasi, G. Geometry optimization in redundant internal coordinates. J. Chem. Phys. 1992, 96, 2856-2860. 40 Baker, J. Techniques for geometry optimization: A comparison of Cartesian and natural internal coordinates. J. Comp. Chem. 1993, 14, 1085-1100. 41 Peng, C.; Schlegel, H. B. Combining synchronous transit and quasi-newton methods to find transition states. Isael J. Chem. 1993, 33, 449-454. 42 Peng, C.; Ayala, P. Y.; Schlegel, H. B.; Frisch, M. J. Using redundant internal coordinates to optimize equilibrium geometries and transition states. J. Comp. Chem. 1996, 17, 49-56. 43 Dewar, M. J. S.; Thiel, W. Ground states of molecules. 38. The MNDO method. Approximations and parameters. J. Am. Chem. Soc., 1977, 99, 4899-4907. 44 Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J. Am. Chem. Soc., 1985, 107, 3902-3909. 45 Stewart, J. J. P. J. Mol. Model., 2007, 13, 1173-1213. 46 Foresman, J. B.; Frisch, AE. Exploring Chemistry with Electronic Structure Methods; Gaussian, Inc., Wallingford, 2015. 47 Birkholz, A.B.; Schlegel, H. B. Using bonding to guide transition state optimization. J. Comp. Chem. 2015, 36, 1157-1166. 48 Austin, A. J.; G. A. Petersson; Frisch, M. J.; Dobek, F. J.; Scalmani, G. Throssell, K. A Density Functional with Spherical Atom Dispersion Terms. J. Chem. Theory and Comput. 2012, 8, 4989-5007.


Page 29 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents Only


Efficient Geometry Minimization and Transition Structure Optimization

Recommend Documents