Efficient Parallel Computation of Molecular Potential Energy Surfaces ...

70

IEEE TRANSACTIONS ON NANOTECHNOLOGY, VOL. 10, NO. 1, JANUARY 2011

Efficient Parallel Computation of Molecular Potential Energy Surfaces for the Study of Light-Induced Transition Dynamics in Multiple Coordinates David Mokrauer, C. T. Kelley, and Alexei Bykhovski

Abstract—The conformational dynamics of molecules that arise due to light-induced transitions are critically important in many biochemical reactions, and therefore dictate the functionality of many types of biological sensors. Therefore, researchers of biological science and biological-inspired technology often need to prescribe the molecular geometry of the stable states and the associated transition trajectories that occur as a result of external excitation, e.g., light-induced transitions from the ground state to the excited states. The traditional approach to study this type of phenomenology is to limit the number of varying molecular coordinates to one or a few due to the considerable computational expense of the required physically modeling required for generating an accurate physical model for analysis. While the conformational dynamics for some types of simple molecules (e.g., retinal) are known to be adequately described by one or few numbers of molecular coordinates, light-induced transitions in arbitrarily complex molecules can be expected to involve the influence of multiple coordinates, and their influence can be expected to vary as a function of time. The research reported here will address the development of parallel computational algorithms that allow for the highly efficient study of molecular conformational dynamics over multiple numbers of multidimensional energy surfaces. Here, the goal is the development of a simulation tool that is capable of: constructing physically accurate multidimensional potential energy surfaces (i.e., from first-principle physical modeling codes); deriving the natural trajectories to local minima within individual surfaces; and that allows for dynamics human interfacing for specifying the transition between energy surfaces and the number of coordinates to be used for the optimization within a particular energy surface. As will be illustrated, this type of physics-based simulation tool will allow researchers to efficiently explore the light-induced conformation dynamics associated with complex biomolecules, and therefore, be a useful tool for the design of biological-sensing processes in the future. Index Terms—Algorithm design and analysis, chemical sensors, 1H NMR, nanoelectonics, stibenes, UV-vis spectroscopy.

Manuscript received February 9, 2010; accepted June 17, 2010. Date of publication July 15, 2010; date of current version January 26, 2011. This was supported was by US Army Research Office under Grant W911NF-07-1-0112. The review of this paper was arranged by Associate Editor D. L. Woolard. D. Mokrauer and C. T. Kelley are with the Department of Mathematics, North Carolina State University, Raleigh, NC 27695 USA (e-mail: [email protected]; [email protected]). A. Bykovski is with the Department of Electrical and Computer Engineering, North Carolina State University, Raleigh, NC 27695 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TNANO.2010.2058862

I. INTRODUCTION OLECULES can transit from one stable (or metastable) state to another due to the excitation caused by the introduction of light [4], [5], [8]–[10]. Understanding this natural phenomenon would help a broad range of new sensors to be built. The occurrence of these transitions in a managed environment can be used to indicate the presence of other substances. The geometry of the molecule can trigger a chain of biochemical reactions in its surroundings, thus, knowing the geometry is imperative. In order to completely understand these transformations, it is useful to know not only the geometry of the initial stable state, but the geometry of the final state as well. These states exist as minima on a potential energy surface (PES). To fully comprehend this transition, it is necessary to possess the path on the PES that connects these two minima. The number of coordinates it takes to fully describe a molecule’s geometry is 3N−6, where N is the number of atoms. Considering that the PES is a function of those 3N−6 coordinates, it would be virtually impossible to compute all the points on the surface for even small molecules. A hypercube with k gridpoints per dimension would require k 3N −6 evaluations. One method often used by researchers to simplify this complexity is to constrain the problem to only a single dimension or perhaps 2-D [6]. If the practitioner has chosen the dominant coordinates, then a good approximation of the transition path and barrier may be obtained. Searching in only 1-D or 2-D assumes that all other coordinates remain nearly constant throughout the entire process, and unfortunately it may be difficult to identify, which coordinates are the most crucial if this is a proper assumption. Since exploring all dimensions is infeasible, we aim to identify the dominant coordinates throughout the transformation. The set of dominant coordinates at one moment may not be the same as at a different moment. Software that allows a user to interface during a study of a transformation on a specific subset of coordinates so that they may dynamically change the search coordinates at opportune moments in the transformation, would give a much more accurate approximation of the path. The user’s experience with simple molecules along with their knowledge of the geometry of the molecule at specific points in the transition would give them an advantage in anticipating, which coordinates may be the dominant set during any section of the path. This initial study seeks to identify algorithms that efficiently take advantage of high performance computing infrastructure to compute portions (or all) of the PES for a specific set of coordinates. Even with the simplification of the problem to a small set

M

1536-125X/$26.00 © 2010 IEEE

MOKRAUER et al.: EFFICIENT PARALLEL COMPUTATION OF MOLECULAR POTENTIAL ENERGY SURFACES

of coordinates, there is a significant computational burden. If we would like to evaluate k gridpoints in n dimensions that requires k n computations, the potential energy at each of these points is computed using computational chemistry software. We will show that these points cannot be computed trivially in parallel on our initial test molecule stilbene. Since these computations for a set of selected coordinates are constrained optimizations on a grid, we cannot simply edit the value of the chosen coordinates and submit. We exercised this algorithm on a 2-D problem that exhibited the characteristics we would like. Stilbene has the potential to be used in different biomolecular-sensing applications, making it a good candidate with which to begin development. The algorithm that we developed does successfully compute an entire PES for stilbene in 2-D. Once a full PES has been obtained in 2-D, we seek to compute a path between local minima on the ground state by simulating light-induced excitations. Since we would like to simulate the natural path the coordinates would follow to minima, we chose the optimization method to be continuous steepest descent. [1] Using the stilbene data we calculated with our new algorithm, we attempted to simulate between two known stable states. In the 2-D that we chose, no feasible path existed. A molecule with a known transition path between two stable states is 2butene [6]. We applied our algorithm and drew a 2-D PES for 2-butene, and then were able to successfully approximate the 2-D transition between minima. II. MOTIVATION The PES associated with the ground state and numerous excited states were obtained using Gaussian 03 [7]. The electron transition information was obtained using the Fermi golden rule [2] in the time-dependent Hartree–Fock (TDHF) and timedependent density functional theory (TDDFT) approximations. Case studies performed in this paper involve ground and several lowest excited energy surfaces of single molecules 2-butene and stilbene. Both molecules are known to have tran- and cisisomers and can interconvert under the influence of light. The stilbene molecule is of particular interest since it can be used in various applications and fields, particularly in biological sensors. For example, different functionalized forms of stilbene (stilbene derivatives) were used for conjugation with DNA [13] and proteins [11] in sensing studies. In the case of the protein work, stilbene derivatives were used to demonstrate a bluefluorescent optical sensor. The main isomerization coordinate of both stilbene and 2-butene is a torsion over the C=C double bond. Torsion angles are determined by four atoms and the three bonds connecting them. A torsion angle is varied by rotating the middle bond of the 3. Fig. 1(a) illustrates the coordinate as the angle between the atoms labeled 1, 2, 3, and 4, which will vary by rotating the bond between 2 and 3. This coordinate will be denoted as D1 in later plots containing potential energies of stilbene. Both molecules are relatively rigid: rotation about the C=C bond is a high-energy process (on the order of 2 eV). The second coordinate of stilbene considered in this paper describes a symmetric rotation of phenyl groups on both carbon atoms of the double bond (labeled D2 in the graphs), which requires

71

Fig. 1. Stilbene torsion angles describing possible transition between stable states(cis to trans). (a) Cis-stilbene with main torsion angle labeled D1. (b) Rotation of D1 toward trans-stilbene. (c) Secondary torsion angle defining phenyl group labeled D2. (d) Rotation of D2 to trans-stilbene.

lower energy. D2 is illustrated in Fig. 1(c) as the torsion angle between atoms 1, 2, 3, and 4. In this paper, we studied the energies of stilbene and 2-butene with the hybrid B3LYP (TDDFT) functional and split-valence electron basis with effective core potentials [12] and ab initio with Hartree–Fock. The 2-butene’s total energy diagram as a function of its main isomerization coordinate (C=C double bond) was studied previously ab initio using Hartree–Fock, and the switching behavior was demonstrated for the excitation of 0–1 electronic transitions by UV light close to 8.5 eV in tran- and cis- 2-butene [6]. In the present study, these earlier results were confirmed by the TDDFT simulations with the hybrid B3LYP functional and splitvalence electron basis with effective core potentials [12]. The excitation energy in the new studies were found to be 7.29 eV for trans- 2-butene, and 7.35 eV for cis- 2-butene, which are marginally less than the earlier work. Moreover, the 2-butene energy was studied with TDDFT in the space of the two internal coordinates. In addition to the isomerization coordinate (C=C double bond torsion), a second coordinate considered is the rotation of tail CH3 groups (i.e., D2 coordinate clarified in Fig. 2). III. INITIAL ATTEMPT AT PARALLELISM We begin with a single stable geometry and select a pair of coordinates. Each point on the 2-D grid is evaluated separately by simply perturbing the two coordinates from this initial geometry. In this paper, steps between gridpoints were 10◦ and the coordinates we chose are torsion angles, which have a range of 360◦ , therefore, there are a total of 372 = 1369 separate computations. These 1369 computations are submitted separately to a high-performance computing cluster, and results were generated in a reasonable amount of time (less than 1 day). Fig. 3 is a scatter plot of the result of this approach for the ground state of stilbene, where every point was generated from a single conformation. The data have some very large jump discontinuities. These are infeasible since the surface should be continuously

72


Fig. 2. 2-butene torsion angles describing possible transition between stable states(cis to trans) (a) Cis-2-butene with main torsion angle labeled D5. (b) Rotation of D5 toward trans-2-butene. (c) Secondary torsion angle defining phenyl group, labeled D8. (d) Rotation of D8 to trans-2-butene.

Fig. 3. Stilbene ground-state PES generated by editing coordinates within a single geometry.

differentiable. The picture also seems to have some small clusters of outliers, which are less of a problem because we can interpolate without including them. The problem with this approach is that the energy values on the PES are computed with a constrained optimization, which does not yield accurate results for mapping the PES. More specifically, the internal optimizations within the computational chemistry software are trust region algorithms, however, while any quasi-Newton method will converge to a minimum with a good initial iterate [3], this approach does not use a good initial iterate for all of these points. By choosing the same initial geometry and varying only two coordinates, we have nearly used the same initial iterate for every point on the surface. If this was a good assumption, then the geometry of the molecule would not be changing much throughout the surface, and thus, the coordinates we have selected are of little importance. The solution to this problem is to use a form of continuation. We will show that a feasible PES can be generated by ensuring that the points have good initial iterates that lead to accurate results. IV. EXPANDING PERIMETER We used a succession of gridpoints in order to simultaneously manage good iterates and take advantage of parallelism. Each

Fig. 4.

2-D expanding perimeter file dependence.

point on the grid, assuming it converges successfully, should be a good initial iterate for any points surrounding it. This is due to the fact that the surfaces are believed to be continuously differentiable. The scheme is illustrated in Fig. 4. Since the user begins with one geometry (in our case it is a local minimum), that conformation serves as the initial iterate for the eight points in a rectangle surrounding it on a 2-D grid (labeled 0 through 8). Once that round of optimizations have all returned (converged or not), a new rectangle is produced. Each optimized geometry in the finished round will serve as the initial iterate for a new point on the next rectangle. The corners of the rectangle are the initial iterates for three points on the next round (e.g., point 1 is the iterate for points 9, 10, and 24). In the event that a file does not converge, the nearest point of a lesser number replaces it as the initial iterate for the point (e.g., if point 12 did not converge, then point 11 will be the initial iterate for both 28 and 29). Each round has eight points more than the previous round. Each point in a round is submitted independently as a Gaussian job, and this parallelism maintains a continuation that leads to accurate results. In order to successfully evaluate the energy near local maxima on the PES, we needed to also utilize some of the symmetry we expect in the final geometry. Atoms 2 and 3 in Fig. 1(a) are bonded to both a ring and a single hydrogen atom. By preprocessing the initial iterate so that the hydrogen atom lay 180◦ opposite the ring, we were able to draw the full PES. Fig. 5 is a scatter plot generated by the expanding perimeter algorithm with preprocessing. There are a few outliers in the scatter plot, which I removed by visual inspection. Fig. 6 is the resulting cubic spline interpolation of the scatter plot with those outliers removed, and Fig. 7 includes a few excited states above the ground state. V. OPTIMIZATION Our goal is to find a path from a local minimum on the ground state to a different local minimum on the ground state. The optimizer will excite the molecule to an excited state, and then minimize on that excited state. Next, we minimize from that

MOKRAUER et al.: EFFICIENT PARALLEL COMPUTATION OF MOLECULAR POTENTIAL ENERGY SURFACES

Fig. 5. Scatter plot of Stilbene ground state generated via expanding perimeter from a single point.

Fig. 6. Stilbene-ground-state interpolated PES generated via expanding perimeter.

Fig. 8.

73

Unsuccessful 2-D optimization of stilbene.

tion solver with ODE: u˙ = − f (u). In this case, the gradient is a finite difference computed from the interpolation. This optimizer is a natural choice since we wish to mimic the natural trajectory. The ODE solver, in this case ODE45 in MATLAB, will manage the error and stepsize, thus, it is unlikely to cross a barrier to find a minimum. Using the gridpoints computed by the expanding perimeter algorithm, we interpolate using a cubic spline. The explicit algorithm looks as follows: 1) input x0 ; 2) select state0 ; 3) x1 = minimum on state0 from initial value x0 ; 4) while x0 − minimum on ground state from point xi ≤ & i < maxit a) i = i + 1 b) Select statei c) xi = minimum on statei from initial value xi−1 . We initially applied the optimization algorithm to the surfaces generated for stilbene. Unfortunately, no matter which state we chose for the first excitation, there was never a path leading away from the original minimum on the ground state. Fig. 8 is a sample of an optimization of stilbene, which returns to the exact location from which it began. This led us to find a different molecule that would verify our concept. 2-butene is a molecule with a known path between local minima via lightinduced excitation [6]. Fig. 9 displays a successful optimization, where a path has been generated for 2-butene.

VI. HIGHER DIMENSIONS Fig. 7.

Stilbene ground state along with excited states 3, 5, and 7 interpolated.

point on the ground state to see if it relaxes to a different point. If not, we choose a different excited state and minimize again. The user is interfacing with the optimizer and selecting the states that are transitioned to. Each time a local minimum is found on that state, the user selects a new state. The minimizations on each state are performed using continuous steepest descent [1]. Continuous steepest descent uses an ordinary differential equa-

Up to now, all of the results we have provided have been verifications of known results. In Fig. 10, we have an example of a 1-D PES with two local minima. The biochemical compound in the figure consists of a DNA fragment chemically bonded to trimethoxystilbene (TMS) carboxamide. Although, we know that two stable states exist researchers have been unable to describe the transformation between these two states. We expect that by properly exploring a large set of coordinates, we will be able to correctly simulate this transition.

74


VIII. CONCLUSION We have developed an algorithm that allows a user to compute the optimal geometry and energy for nearly any feasible point on the PES using ab initio computational chemistry software. Using this algorithm, we have also shown that paths exist that transform molecules from one stable conformation to another in 2-D through excitation. The algorithm extends indefinitely into multidimensional space. This software will allow researchers to predict and harness light-induced transformations, which will lead to new sensing technology. D. Mokrauer would like to thank Dr. D. L. Woolard, of the US Army Research Office, for the mentoring assistance in the preparation of this paper. REFERENCES Fig. 9.

Successful 2-D optimization of 2-butene.

Fig. 10. Predicted energy of the double-strand DNA fragment TGCGCA. One DNA strand is capped with TMS. The positions of stable conformations are highlighted with arrows. The zeroth, first, second, and ninth states are shown.

VII. RESOURCES All computations were performed on the high-performance computing cluster at North Carolina State University. Our chassis has 56 dual core Xeon processors with 2 GB distributed memory per core and dual gigabit ethernet interconnects. The operating system is Red Hat Linux 2.6.9. Potential energy computations are performed using Gaussian 03. Script editing is done with both Gawk 3.1.3 and Python 2.5.4. Figures and optimizations were produced in MATLAB 7.8.0.347.

[1] R. Courant, “Variational methods for the solution of problems of equilibrium and vibration,” Bull. Amer. Math. Soc., vol. 49, pp. 1–43, 1943. [2] S. Gasiorowicz, Quantum Physics. New York: Wiley, 1974. [3] C. T. Kelley, Iterative Methods for Optimization. Philadelphia, PA: SIAM, 1999. ´ Kvaran, B. I. Asgeirsson, ´ [4] A. and J. K. F. Geirsson, “1H NMR and UV-vis spectroscopy of fluorine and chlorine substituted stilbenes: Conformational studies,” J. Mol. Struct. vol. 563/564, pp. 513–516, 2001. [5] F. D. Lewis, “Formation and reactions of stilbene exciplexes,” Acc. Chem. Res., vol. 12, pp. 152–158, 1979. [6] Y. Luo, B. L. Gelmont, and D. L. Woolard, “Bio-molecular devices for terahertz frequency sensing,” in Molecular and Nano Electronics: Analysis, Design and Simulation, J. M. Seminario, Ed. ch. 2, Elsevier, 2007, pp. 55–81. [7] F. R. Clemente, M. J. Frisch, A. E. Frisch, and G. W. Trucks, Gaussian 09 User’s Reference. Wallingford, CT: Gaussian, Inc., 2009. [8] I. J. Palmer, I. N. Ragazos, F. Bemardi, M. Olivucci, and M. A. Robb, “An MC-SCF study of the S1 and S2 photochemical reactions of benzene,” J. Amer. Chem. Soc., vol. 115, pp. 673–682, 1993. [9] J. Quenneville and T. J. Martinez, “Ab initio study of cis-trans photoisomerization in stilbene and ethylene,” J. Chem. Phys., vol. 107, pp. 829– 837, 2003. [10] J. Saltiel, A. S. Waller, and D. F. Sears, Jr., “The temperature and medium dependencies of cis-stilbene fluorescence: The energetics for twisting in the lowest excited singlet state,” J. Amer. Chem. Soc., vol. 115, pp. 2453– 2465, 1993. [11] A. Simeonov, M. Matsushita, E. A. Juban, E. H. Z. Thompson, T. Z. Hoffman, A. E. Beuscher IV, M. J. Taylor, P. Wirsching, W. Rettig, J. K. McCusker, R. C. Stevens, D. P. Millar, P. G. Schultz, R. A. Lerner, and K. D. Janda, “Blue-fluorescent antibodies,” Science, vol. 290, pp. 307– 313, 2000. [12] W. J. Stevens, H. Basch, and M. Krauss, “Compact effective potentials and efficient shared-exponent basis sets for the first and second-row atoms,” J. Chem. Phys., vol. 81, pp. 6026–6033, 1984. [13] J. Tuma, R. Paulini, J. A. Rojas Stutz, and C. Richert, “How much pistacking do dna termini seek? Solution structure of self-complementary dna hexamer with trimethoxystilbenes capping the terminal base pairs,” Biochemistry, vol. 43, pp. 15680–15687, 2004.

Authors’ photographs and biographies not available at the time of publication.