An Efficient Strategy to Generate Atom Connecting Positions

TL;DR

This paper introduces a fast, automated method for generating atom connecting positions (ACPs) using spherical optimization and VSEPR theory, improving efficiency and accuracy in transition state and minimum energy path calculations.

Contribution

The paper presents a novel, automated approach for ACP generation that eliminates manual effort and enhances robustness in computational chemistry workflows.

Findings

Method is more efficient than manual generation.

Approach is robust across multiple examples.

Automates ACP generation without manual interference.

Abstract

Atom connecting positions(ACPs) are positions where an atom is connecting to another one or a few atoms, which is needed when constructing final state used in chain-of-state(CoS) methods for transition state(TS) locating and minimum energy path(MEP) searching, especially with bond formation. However, ACPs are generated with chemical insight and experience, which is not only low efficient and time wasting, but the manually generated structure may be far from the optimized one. A efficient method is presented here for generating ACPs which is based on spherical optimization and VSEPR theory without manual interfering. Several examples are testified to prove the efficiency and robustness of the method.

Tables

Table 1. Table 1 : Angle Difference between VSEPR-SOPT and DFT Optimization

Item	Angle Difference
\ceAu20@Aua	${0.4}^{\circ }$
\ceAu20@Aue	${13.43}^{\circ }$
\ceAu20@Auc	${2.24}^{\circ }$
CNT@circle	${2.48}^{\circ }$
\ceC60@circle	${8.13}^{\circ }$
\ceAu1/\ceCeO2	${13.19}^{\circ }$

Equations

$$\begin{array}{rl}& \underset{\vec{x}}{\min }E=E(x_1,x_2,...,x_n)\\ & \\ & s.t.(x_1-y_1{)}^2+(x_2-y_2{)}^2+...+(x_n-y_n{)}^2=R^2\end{array}$$

$$E^{\prime }=E(x_1,x_2,...,x_{n-1},f(x_1,x_2,...,x_{n-1},R))$$

$$\begin{array}{rl}{F}_i^{\prime }=-\frac{\partial E^{\prime }}{\partial x_i}=& -(\frac{\partial E}{\partial x_i}+\frac{\partial E}{\partial x_n}\frac{\partial x_n}{\partial x_i})\\ =& F_i-F_n\frac{x_i-y_i}{x_n-y_n}\end{array}$$

$$\begin{array}{r}E(\vec{x})=\sum _{l=1}^N\frac{1}{|\overrightarrow{r_l}-\vec{x}{|}^4}\end{array}$$

$$\begin{array}{r}E(\vec{x})=\sum _{i,j,k=-1}^1\sum _{l=1}^N\frac{1}{|\overrightarrow{r_l}+i\vec{a}+j\vec{b}+k\vec{c}-\vec{x}{|}^4}\end{array}$$

$$\begin{array}{r}E(\overrightarrow{x_1},\overrightarrow{x_2},...,\overrightarrow{x_M})=\sum _{i=1}^N\sum _{j=1}^M\frac{1}{|\overrightarrow{r_i}-\overrightarrow{x_j}{|}^4}+\sum _{i,j=1,i\ne j}^M\frac{1}{|\overrightarrow{x_i}-\overrightarrow{x_j}{|}^4}\end{array}$$

Full text

DOC, Tsinghua University]Department of Chemistry, Tsinghua University, Beijing, 100084, P.R. China

DOC, Tsinghua University]Department of Chemistry, Tsinghua University, Beijing, 100084, P.R. China

DOC, Tsinghua University]Department of Chemistry, Tsinghua University, Beijing, 100084, P.R. China

1 INTRODUCTION

Initial and final state should be prepared before doing transition state(TS) locating and minimum energy path(MEP) searching using chain-of-state(CoS) methods like nudge elastic band(NEB) method 1, 2, 3, 4, 5, 6. Generally speaking, with a bond formation described in the MEP, the initial state is obtainable, whereas the final state is usually unknown and should be constructed artificially. Since there is bond reforming and atom transferring, determining suitable positions for the atom to be positioned is really important. These positions are defined as atom connecting positions(ACPs). However, these positions are generally determined manually based on chemical insight and experience, because these positions correlate with geometry structure as well as electronic structure.

Valence shell electron pair repulsion(VSEPR) theory is a good start, since the ACPs determined by chemical insight are generally based on the repulsion of connected atoms. When one atom is connecting to the destination atom, suiting the destination atom to one of the VSEPR shapes like Linear, Tetrahedral, Pentagonal planar or Octahedral is a possible option. ACPs could then be generated with certain VSEPR shape.

But direct usage of VSEPR shape model is not only limited and complicated, but also would fail under some situations. Firstly, there are too many VSEPR shape models and variations7. especially when distortion exists. Next, VSEPR shape model can only be used when the transferred atom connect to another atom, but it could also connect with multiple ones simultaneously. Besides, VSEPR would fail when the atoms are not connected by chemical bond, like metal surface and cluster, in which atoms are packing together. Choosing connected atoms for suiting the shape is very important for VSEPR, but it’s difficult to classify two atoms to be connected or not when they are weakly connected, and different choice may give out different result. All these drawbacks may downgrade the robustness of the method. Therefore, we need to purpose a general, efficient and robust method for determining ACPs connecting with not only specific atom as a center, but a virtual center like the mid-point of two atoms or the center of a benzene ring. The center is defined as kernel.

By rethinking the motivation of VSEPR theory, we found that the main idea is to find out the low-repulsion positions around the spherical surface centered at kernel. Therefore, we believe the ACPs are actually the collection of low-repulsion positions. This leads us to a constrained optimization problem. By optimizing all the positions on the spherical surface centered at the kernel, we may get all ACPs. By introducing spherical optimization(SOPT) method 8, we formulated the framework of VSEPR-SOPT, which turns the chemical insight into a robust and efficient algorithm.

2 METHODS

2.1 Spherical Optimization

The detail of spherical optimization(SOPT) has been described by Abashkin, Y. and Russo, N. 8, so it will be introduced briefly. Considering two atom system with $N$ atoms and noted as $x$ and $y$ with root-mean-square deviation of distance as $R$, the constrained optimization problem is described as

[TABLE]

where $n=3N$ is the total degrees of freedom. $x$ and $y$ are also named as target and anchor, since $x$ will be optimized and $y$ is fixed during the optimization. The energy function is rewritten by including the constraint and eliminate $x_{n}$,

[TABLE]

and the force is rewritten as

[TABLE]

where $F$ is the force obtained from electronic structure calculation, and $F^{\prime}$ is a $n-1$ with $x_{n}$ eliminated. Thus, we convert this particular constrained optimization with $n$ variables to a regular optimization problem with $n-1$ variables. More details of the algorithm are neglected and can be found elsewhere.

2.2 VSEPR-SOPT Model

Combining SOPT with VSEPR theory is tricky, since the target and anchor system are not real system, and the energy function remains unknown.

Firstly, let’s consider locating one ACP on the sphere centered at kernel. Locating multiple ones can be regarded as locating ACP one after another. Positions on the sphere centered at kernel with radius $R=L$ is the collection of possible positions. Here $L$ is the length between kernel and the position on the sphere, which is given by user or get from database. So the target system is constructed with a pseudo atom positioned on the sphere, and the anchor system contains only the position of kernel.

The energy of the system is to describe repulsion between the pseudo atom and any other atoms in the target system. Minimizing the energy will give out the ACP directly. But the formula of energy is hard to determine since the system is not real, and VSEPR repulsion has not been formulated. Besides, we expect the calculation to be as fast as possible without introducing too much parameters. Inspired by Lennard-Jones potential, we believe the repulsion energy between a pair of atoms can be formulated as $\frac{1}{|\vec{r}|^{n}}$, where $n$ is a constant number. It turned out $n=4$ will give out the best results. For molecular system, the formula of energy is

[TABLE]

where $N$ is the total number of atoms in real system, $\vec{x}$ is the position of pseudo atom and $\vec{r_{l}}$ is the position of $l^{th}$ atom in real system. It’s worth mentioning that no external parameter is introduced in the function. It seems weird, but the results turned out to be good enough.

For periodic system, the contribution of replicas should be counted as well. Therefore, another 8 neighbor cells are introduced. The energy function is

[TABLE]

where $\vec{a},\vec{b},\vec{c}$ are unit vectors of the cell.

For $M$ ACPs, we just need to construct the target system contains $M$ pseudo atoms, and duplicate kernel $M$ times in anchor system, adjust the radius of the sphere $R=L\sqrt{M}$. Interactions between pseudo atoms should be added to the energy function

[TABLE]

The rest part remains the same. And for periodic system with multiple ACPs, the formula will be modified as presented above. With the formula of energy function, SOPT now is able to be utilized for finding ACPs. Since the simplicity of the procedure, the execution of the algorithm is really fast.

2.3 Position Sampling on Spherical Surface

In order to acquire all ACPs on the spherical surface, we need to sample all positions on the surface. Therefore, a sampling algorithm should be introduced. Using polar coordinate system and sampling $\theta$ and $\phi$ uniformly is simple and convenient, but the sampling density around the pole is much higher than that around the equator. However, a uniform sampling on the spherical surface method is expected. Here we use a method introduced by Cory Simon9. By generating three standard normally distributed numbers $X$, $Y$, $Z$ to form a vector $V=[X,Y,Z]$ and normalize it, the vector is uniformly distributed on the surface of sphere. Normally, sampling 100 times is adequate for covering the entire surface, but the number can be increased if necessary.

2.4 Flowchart of VSEPR-SOPT

We provide the flowchart for the algorithm of VSEPR-SOPT below,

Input the system, kernel position, number of ACPs needed($N_{n}$). 2. 2.

Add a pseudo atom to the system and construct energy function. 3. 3.

Sampling on the sphere $N_{s}$ times with the algorithm described above, optimize the pseudo atoms system with the potential described in Eq.6 with SOPT. 4. 4.

Repeat step 2-3 $N_{n}$ times to get all ACPs needed.

3 RESULTS

Here we use the notation A@B to represent system A with B as the kernel, e.g. \ceCH2O@C means a molecular system \ceCH2O, with the carbon atom as the kernel. It’s worth mentioning that the center of a ring can be kernel as well. The kernel will be noted as circle in these situations. In the following, we tested several kinds of system, including small molecules (\ceCH4, \ceCH2O), clusters (\ceAu20, \ceC60), complicated molecular system(complicated \ceAu6 cluster) and heterogeneous system(\ceCeO2-based single atom catalyst(SAC)). In the related figures, ACPs are labeled with pseudo atoms in violet, and only kernel(or kernel related) atoms and pseudo atoms are in ball-stick format and other atoms are in line format. The B3LYP functional with default parameters in Gaussian 0910 are used in all DFT calculations and geometry optimization if not specified. 6-31G(d) basis sets are used for C, O, H and LANL2DZ are used for Au, Fe. Difference of angle between VSEPR-SOPT and DFT geometry optimization are calculated with other atoms aligned, and the numbers are shown in Table.LABEL:table:angle-diff. This proves the validity of our method.

3.1 \ceCH4 and \ceCH2O

System

VSEPR-SOPT results of \ceCH4@C, \ceCH4@H are shown in Fig.1(a) and 1(b). For \ceCH4@C, all 4 face centers of the tetrahedral are acquired, and for \ceCH4@H, linear shape is implied and the other end is acquired as ACP. Results of \ceCH2O@C, \ceCH2O@O are shown in Fig.2(a) and 2(b). For \ceCH2O@C, a trigonal bipyramidal is displayed and the positions off the \ceCH2 plane are acquired. However, for \ceCH2O@O, the lone-pairs’ positions acquired violates the chemical rule, since we expect the lone-pairs to be in the plane of \ceCH2O, but the positions given by VSEPR-SOPT are off the plane. This is due to the free sampling of the spherical surface and the absence of electronic structure. But the result is still acceptable. The results shown here satisfy chemical insights in general, and the weakness is endurable.

3.2 \ceAu20 Cluster System

\ce

Au20 cluster11 has 3 kinds of Au: Au along the edges(Aue), Au at the apexes(Aua) and Au at the center of each face(Auc). ACPs for \ceAu20 with each kind of Au atom acquired by VSEPR-SOPT are shown in Fig.3(a)-3(c). For each kind of Au atom, only one position is identified as ACP. To show the correctness of the result, we further testified geometry structure of CO absorbed on \ceAu20 with DFT calculation, as shown in Fig.4(a)-4(c). The results given by DFT match the VSEPR-SOPT results astonishingly. The angle differences are $0.4^{\circ}$, $13.43^{\circ}$, $2.24^{\circ}$ for \ceAu20@Aua, \ceAu20@Aue, \ceAu20@Auc, respectively. Relative large angle difference for \ceAu20@Aue is account for distortion of \ceAu20 cluster, but it is still acceptable.

However, suiting these kinds of Au atom with regular VSEPR shape model is hard or even impossible, since these Au atoms are not connecting with chemical bonds but packing together. It’s really difficult to determining the VSEPR shape for these Au atoms. For a Aua atom, 3 Au atoms are coordinated, while the angles are only $60^{\circ}$ instead of $120^{\circ}$ to suit a Triangle model or $90^{\circ}$ to suit a Octahedral model. For a Aue atom, 6 Au atoms are coordinated in a strange shape which doesn’t belong to any VSEPR model. And Auc atom is coordinated with 9 Au atoms, which is far beyond the range of VSEPR shape model. This may cause failure for regular VSEPR model, whereas the result given by VSEPR-SOPT match our chemical insight and the reality remarkably.

3.3 \ceC60 and Carbon Nanotube(CNT)

System

For \ceC60 and CNT, all carbon atoms are equivalent, but the ACPs can be inside or outside the cage/tube. For simplicity, only outside of cage/tube is testified. In these cases, center of a carbon circle can also be kernel, and these circles are equivalent as well. So we use both carbon and the center of circle to be the kernel.

\ce

C60 has two kinds of carbon circle, \ceC5 circle and \ceC6 circle. Here we only use \ceC6 circle and \ceC6 could be processed with the same strategy. Using VSEPR-SOPT, we can identify the ACPs on \ceC60 easily, as shown in Fig.5(a) and 5(b). \ceTiC60 presented by Sun, Qiang, et al. 12 is optimized with PAW method implemented with VASP. The result shown in Fig.5(c) agrees with the \ceC6 kernel situation very well and the angle difference is about $8^{\circ}$.

As for CNT, the same setups are executed. The VSEPR-SOPT results with carbon and \ceC6 as kernel are shown in Fig.6(a) and 6(b), Fe-CNT system is optimized with DFT and presented in Fig.6(c). DFT result shows that Fe is just above the center of the \ceC6 circle, as what we acquired with VSEPR-SOPT. The angle difference is less than $3^{\circ}$, which agrees with the reality extraordinarily.

For \ceC60 system and CNT system, ACPs acquired by VSEPR-SOPT are not only chemically meaningful, but also satisfy the DFT calculation results. Whereas for regular VSEPR shape model suiting, it would be really hard since the kernel here is not an individual atom but the center of a carbon circle.

3.4 A Complicated \ceAu6

Cluster

\ce

Au6 cluster13 is a very large cluster system, which contains 6 equivalent Au atoms, 6 P atoms, 2 N atoms and 20 phenyl groups(Fig.7). Environment of the Au atom is very complicated. One Au atom is not only connecting with one N and P, but connecting with other Au atoms weakly, the phenyl groups around the atom have influence as well. Therefore, it’s impossible to determine the VSEPR shape of Au atom. It’s still very hard to get the positions manually with chemical insight. However, with the help of VSEPR-SOPT, ACPs could be acquired easily, ACPs acquired by the algorithm are shown in Fig.8. Noted that the VSEPR-SOPT does not need the connection relationship at all, which is critical for regular VSEPR shape model. This example proves the robustness of the method.

3.5 \ceCeO2-based Single Atom

Catalyst

\ce

CeO2-based single atom catalyst could be used as heterogeneous catalysis, like \ceCO2 reduction and CO oxidation. \ceAu1/\ceCeO2 has been reported14 and the geometry optimizations are implemented with VASP, as shown in Fig.9. CO can be regarded as ACP detector in this example. It’s worth mentioning that this system is a periodic system. So the energy formula of periodic system should be used. The VSEPR-SOPT result is shown in Fig.10. And the ACP acquired agrees well with the position of carbon mentioned in the article. Angle difference is $13.19^{\circ}$, which is a really good result, Showing the validity of the method in heterogeneous systems.

4 SUMMARY AND DISCUSSION

In this article, we purposed an efficient method for locating atom connecting positions(ACPs). The basic idea is to substitute chemical insight from VSEPR theory to a constrained geometry optimization problem. With the help of spherical optimization(SOPT), the chosen energy function and sampling method, all ACPs are acquirable and satisfy chemical insight basically. Several typical systems are testified, including small molecules, simple clusters, complicated molecular system and heterogeneous system. In all systems except \ceCH2O@O, chemical insight is satisfied rigorously, while for \ceCH2O@O, chemical insight is partially satisfied as well, but the result is still acceptable. The angle differences are less than $20^{\circ}$, which is astonishingly small. These examples demonstrate the validity and robustness of the method. And the method can be useful in reaction searching and reaction network construction.

{acknowledgement}

This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 21590792, 91426302, and 2143005) to J.L. and the Thousand Talents Plan for Young Scholars to H.X. The calculations were performed using the supercomputers at the Computational Chemistry Laboratory of Department of Chemistry under Tsinghua Xuetang Talents Program. We thank Biao Yang from Tsinghua University for providing the structure of \ceAu6 cluster.