Research: Overview

Research Interests


Modern theoretical methods have advanced to an extent that the microscopic details of chemically important processes can be now be investigated using novel algorithms and high speed computer. My research aims to further develop and exploit these emerging capabilities by devising new molecular dynamics and electronic structure techniques and applying these in chemically and biologically important systems. Specific current projects going on in my group are:

Ab initio molecular dynamics investigations of the addition of organic molecules to semiconductor surfaces.

Ab initio molecular dynamics and path integral studies of the solvation and transport of topological charge defects in hydrogen-bonded liquids including water and aqueous solutions and at the corresponding vapor/liquid interfaces.

Development of new, simple real space basis set approaches for density functional based ab initio molecular dynamics.

Development of novel adiabatic dynamics and variable transformation based methods for conformational sampling, biomolecular structure prediction, and crystal polymorphism exploration.

Development of rigorous theory of molecular pseudopotentials.

Development of theoretical statistical mechanical analysis techniques for non-Hamiltonian dynamical systems used in molecular dynamics calculations.

Development of a large-scale, parallel, object-oriented simulation packages (PINY_MD and OpenAtom) capable of performing force field based, ab initio and path integral molecular dynamics calculations in a variety of ensembles.


Selected results from these studies are summarized below:



Reactions of organic molecules on semiconductor surfaces.

The chemistry of hybrid structures composed of organic molecules and semiconductor or metal surfaces is opening up exciting new directions if molecular electronics and nanoscale devices. Theoretical studies of the addition of 1,3-butadiene and 1,3-cyclohexadiene, both conjugated dienes, to the Si(100)-2x1 surface obtained a distribution of addition products. These include a [4+2] Diels-Alder adduct with a single surface dimer, a [4+2]-like adduct in which the butadiene bridges two dimers within a row, a [4+2]-like adduct, in which the butadiene bridges two dimers in neighboring rows, and a [2+2]-like adduct bridging two dimers within a row. The proportion of each adduct in the distribution is in agreement with STM measurements. Moreover, a common underlying mechanism capable of rationalizing the formation of each adduct was proposed (JACS 127, 1110 (2005) -- see Figure at right). It was shown that modification of the butadiene by substitution of a fluorine for a hydrogen lowers the free energy barrier for removal of the molecule from the surface, a result that has important implications for surface lithography (PNAS 102, 6654 (2005) -- see Figures below). The figure below shows adducts obtained using 1,3-cyclohexadiene on the Si(100)-2x1 surface.


Animations:
1,3-cyclohexadiene forming [4+2] Diels-alder adduct on Si(100)-2x1
1,3-cyclohexadiene forming [4+2]-like inter-dimer, intra-row adduct on Si(100)-2x1
1,3-cyclohexadiene forming [2+2]-like adduct on Si(100)-2x1
1,3-cyclohexadiene forming [4+2]-like inter-dimer, inter-row adduct on Si(100)-2x1

Because Si(100)-2x1 admits so many different products, it is not possible to created ordered layers on the surface using a molecule like 1,3-cyclohexadiene or 1,3-butadiene. If we switch to 3C-SiC(001)-3x2 silicon-carbide surface, the silicon dimers are further apart, which could reduce the number of accessible products in the distribution. Interestingly, however, this surface still allows a number of products besides the Diels-Alder product, including one that involves a sublayer atom. The animations below shows some of these products.

Animations:
1,3-cyclohexadiene forming [4+2] Diels-alder adduct on 3C-SiC(001)-3x2
1,3-cyclohexadiene undergoing hydrogen abstraction on 3C-SiC(001)-3x2
1,3-cyclohexadiene forming a subsurface adduct on 3C-SiC(001)-3x2

However, consider the SiC(100)-c(2x2) surface, which is characterized by well-separated triple-bonded carbon dimers:


These dimers can also act as dienophiles following the more common Woodward-Hoffman rules, and this strongly restricts the product distribution to the Diels-Alder channel. This can be verified by computing the free energy profile for this and other reaction channels, as shown in the figures below [see, also, Zhang and Tuckerman J. Phys. Chem. Lett. 2, 1814 (2011)]. The average distances between dimers in a row and between neighboring rows are 4.39 and 4.97 angstroms, respectively.




The only reaction channel that shows an initial barrier under about 20 kcal/mol is the [4+2] Diels-Alder adduct, suggesting that an ordered layer could be formed using this particular reaction, which could lead to molecular wires or ordered passivating layers.





Anomalous transport of topological defects in aqueous solutions.

The transport of topological defects in hydrogen-bonded liquids, created by the addition or removel of a proton is a ubiquitous processe that is fundamental in numerous biologically and technologically important systems. In water, hydronium (H3O+) and hydroxide -) ions are known to diffuse at anomalously high rates (compared to the self-diffusion of water). The basic picture of this anomalous transport dates back to an idea of C. J. T. von Grotthuss in 1806, in which he posited that the ions are NOT transported hydrodynamically as intact units but rather that a defect diffuses "structurally" through the hydrogen-bond network via a series of proton transfer reactions. Using the methodology of ab initio molecular dynamics, we have carried out investigations of the underlying molecular mechanism of anomalous transport of hydronium (H3O+) and hydroxide (OH-) in water. Our studied yielded a number of key predictions that were subsequently verified by experiment.

For hydronium, the mechanism begins with a decrease in coordination number of a water in hydronium's first solvation shell from 4 to 3. The lost hydrogen bond is an acceptor hydrogen bond. This change in the coordination shell leaves the water molecule with a coordination pattern that more closely resembles that of hydronium that water, thereby "activating" the hydrogen bond between hydronium and the neighboring water. The proton often begins to "rattle" in the hydrogen bond, meaning that it often passes from one water moiety to the other. An actual charge displacement event seems to occur when the original hydronium starts to acquire a fourth hydrogen bond (an acceptor hydrogen bond), thus giving it a coordination pattern like that of water (Click here for an animation of the process). This causes the proton to transfer. In fact, the coordination reduction in the first solvation shell and coordination increase at the hydronium site often occur in rapid succession, thus leading to a nearly "concerted" mechanism. This is illustrated in the sequence of snapshots at the left.



For hydroxide, it has been thought for over a century that the anomalous mechanism of hydroxide transport could be deduced from that of hydronium by simply reversing all of the hydrogen bond polarities. In fact, the chemistry of hydroxide is not the same as that of hydronium, and hydroxide's transport mechanism is, consequently, rather different. The picture at the left shows a sequence of snapshots from our path-integral simulations, together with the electron localization function (blue isosurface). The electron localization function indicates that the lone pairs of hydroxide are not distinct but become delocalized in a ring-like structure, and it is this structure that leads to the unusual coordination pattern around the hydroxyl oxygen. We call this pattern a "hypercoordinated" pattern because it consists primarily of a fourfold, planar structure that is inactive with respect to proton transfer (panel (a)). This structure is, however, somewhat labile, and it can occasionally change from fourfold to a threefold tetrahedral structure (panel (b)). At the same time, the hydroxide ion can also donate a weak hydrogen bond (panel(c)), and the combination of coordination reduction around the hydroxyl hydrogen and donation of a hydrogen bond leaves hydroxide in a coordination pattern that resembles water. In this state, a proton can be transferred from a first solvation shell water (panel(d)). As the proton is transferred, the ring structure is also "transferred" to the nascent hydroxide, and the process becomes a complete one when the nascent hydroxide acquires a fourth hydrogen bond to complete its solvation shell (panel (f)). These findings for hydroxide have since been confirmed from neutron scattering, X-ray scattering, core-hole spectroscopy and FTIR. For both hydronium and hydroxide transport, we have developed a kinetic model based on population correlation functions that allow the time scales associated with both "rattling" and actual charge displacement events to be calculated from a molecular dynamics trajectory (Click here for Physical Review Letters article).




Our studies of proton solvation and transport in hydrogen-bonded liquids such as water, methanol and methanol-water solutions have elucidated microscopic mechanisms, resolved a number of long-standing controversies, and have uncovered a number of general patterns and principles that appear to govern this process in different liquids. We have also studied the liquid/vapor interface of an HCl solution [H. S. Lee and M. E. Tuckerman, J. Phys. Chem. A 113, 2144 (2009)] and showed, in agreement with other studied, that the hydronium has a clear preference for the interface (see image at right). Calculation of the potential of mean force for transfer of the hydronium from the bulk to the surface and just beyond predict a shallow minimum at the surface of approximately 1.3 kcal/mol. Interestingly, in a 50:50 methanol water solution, in which water tends to form small clusters, the hydronium is preferentially located in the "interfacial" at the edge of the water cluster where it meets methanol. Although the charge defect never penetrates into methanol regions, the simulations [Morrone, Haslinger, Tuckerman J. Phys. Chem. B, 110, 3712 (2006)] predict rattling events between water and water and between water and methanol.


Recently, we studied the temperature dependence of the reorientation of the OH- ion basic solution [see Z. Ma and M. E. Tuckerman Chem. Phys. Lett. 511, 177 (2011)]. The reorientation time has also been measured from femtosecond spectroscopy of the charge-transfer-to-solvent transition by Keiding and coworkers [Chem. Phys. Lett. 481, 9 (2009)], and this figure compares to these measurements. What the AIMD simulations reveal is that as the temperature is lowered, the striking decrease in the reorientation time is closely connected to suppression of the structural diffusion mechanism for hydroxide discussed above.



Real-space basis sets for ab initio molecular dynamics.

Real space approaches in electronic structure have typically employed Gaussian basis sets, which, although very efficient, possess many technical complexities. Plane wave basis sets, on the other hand, although very simple and elegant, often suffer from poor convergence. By employing discrete variable representation (DVR) techniques with plane wave basis sets for long range interactions and our novel screening function methodology, we have developed a simple and rigorous real space based approach for density functional based electronic structure calculations which lead to an order of magnitude gain in efficiency over plane wave methods but retain all of the simplicity of the latter. The images at the left show examples of one-dimensional and two-dimensional DVR functions. The picture shows that these functions are localized at specific points of a quadrature grid. As a consequence, localized orbitals can be represented by a small number of such functions, thereby giving rise to sparse matrices and potentially a linear-scaling algorithm. Moreover, it is possible to carry out simulations with fully converged basis sets with low computational overhead. At present, long-range terms in the energy functional are calculated in reciprocal space, which requires two FFTs (one for the density and one for the Kohn-Sham potential), which scale as NlogN. This basis set approach was used in our recent studies of neat water, and the calculations show a significant improvement in the structural and dynamical properites of water with respect to underconverged basis sets. For a more detailed write up, see our article in J. Phys. Chem. A 110, 5549 (2006), or click here for our previous article in Physical Review B.



Spatial-warping transformations and adiabatic dynamics for enhanced conformational sampling on rough energy landscapes.




The conformational sampling problem is one of the computational grand challenges. If solved, problems such as protein and nucleic acid folding will be significantly impacted. By introducing novel variable transformations into the mathematical expression for the statistical mechanical partition function and introducing adiabatic dynamics schemes, we have developed novel molecular dynamics approaches that lead to very large gains in efficiency in the ability to sample statistically relevant conformations of large chain molecules and to compute free energies and other equilibrium properties. To illustrate the approach in a very simple case, consider a double well potential V(x) with wells at +/-a and a large barrier separating them. In order to predict the equilibrium and thermodynamic properties of such a system, we need to compute the canonical partition function:

If we try compute this partition function using thermostatted molecular dynamics (or Monte Carlo), barrier crossing would be infrequent, decreasing exponentially with barrier height). However, since we are only interested in computing an integral, consider changing from the coordinate x to a new coordinate u according to:

where the Boltzmann factor involves a new potential that we call a reference potential. Substituting this into the partition function yields

If we, therefore, choose the reference potential to be the true potential V(x) in the barrier region and 0 otherwise, then the appearance of the difference in the transformed partition function suggests that there will be no barrier in u space. We, therefore, call these transformations spatial-warping transformations, as they warp configuration space in such a way as to remove unwanted barriers. However, the transformation does not alter the partition function; thus, all thermodynamic and equilibrium properties are retained, and there is no need to unbias configurations generated a posteriori. Indeed, if a molecular dynamics calculation is performed using u instead of x and compared to one using x for several barrier heights, the result is the image at the left, which shows x as a function of time for both simulations and the probability distributions that result (left -- no transformations, right -- transformations). The top two rows show barrier heights of 5kT and 10kT, respectively, and the dramatic slowdown of barrier crossing events in normal molecular dynamics is clear from the two pictures. However, with spatial-warping transformations, there is no slowdown at all of barrier crossing for either barrier height, and the probability distribution P(x) is generated directly from the simulation with no a posteriori reweighting of configurations because the partition function is not fundamentally changed.
In order to apply spatial-warping transformations to polymers and biomolecules, we transform directly on the backbone dihedral angles in order to remove barriers in this part of configuration space. It is possible to develop reference potentials capable of removing "intrinsic" barriers that occur in the torsional part of the potential surface and "unanticipated" barriers that arise on the fly due to nonbonded interactions. Click here to read a more detailed writeup. In order to compare this approach with another popular scheme, parallel tempering, we choose a 50-mer alkane chain treated using the CHARMM22 all-atom force field. Parallel tempering is performed using 10 replicas with temperatures between 300 K and 1000 K. The picture on the right shows comparison (left = thermostatted MD, middle = parallel tempering, right = spatial-warping transformations). The top row shows a plot of dihedral angle number (horizontal) vs. the number of barrier crossing events (vertical). The middle row shows a "Ramachandran" scatter plot of the two central dihedral angles (25 and 26), and the bottom row shows the evolution of the end-to-end distance. It is clear that only the spatial-warping transformation scheme (which uses one replica) is able to visit all 9 basins in the Ramachandran plot and generate and end-to-end evolution that has a well-defined average.

We have also shown that coordinate transformations introduced into the partition function can be used to generate multi-dimensional free-energy surfaces in small sets of collective variables useful for describing a particular process. The original method we developed in 2000 (first presented at a CECAM meeting on free-energy methods) and published in 2002 (click here for full article) requires that a transformation exist that explicitly involves these collective variables. The collective variables are subsequently kept at a high temperature, which allows them to cross any barriers in their part of the configuration space, and they are assigned a high mass so that they evolve slowly during the simulation. Under these conditions, the collective variables are adiabatically decoupled from all other variables, which allows the latter to relax nearly instantaneously to the motion of the collective variables, and it can be shown that the free energy surface is


where Ts is the temperature of the collective variables, and Padb is the adiabatic probability distribution for n collective variables to have values s1,...,sn. This generates the free energy surface at the correct temperature T. The method is known as adiabatic free energy dynamics (AFED). The approach is powerful, conceptually simple, and can be rigorously shown to yield the correct free energy surface. Indeed, we have used it to generate free energy surface of solvated dipeptides [J. Phys. Chem. B 109, 4162 (2005)]. However, the need to employ explicit coordinate transformations is a disadvantage. More recently, we simplified this scheme by introducing s1,...,sn as external variables with conjugate momenta p1,...,pn and a harmonic coupling between the external variables and the actual collective variables. We, therefore, take the basic Hamiltonian to be

Here, the first term is the physical Hamiltonian, and the qs are the collective variables. In this way, AFED calculations can be performed without explicit coordinate transformations (see our recent article, J. Phys. Chem. B 112, 15742 (2008)]. It also allows a wider range of collective variables to be studied, since we do not need to devise a coordinate transformation that includes them explicitly. We employed this approach to map out the free energy surface of the alanine dipeptide and N-acetyl-tryptophan-methylamide, as well as an alanine hexamer, shown below. In particular, we focus on the number of backbone hydrogen bonds and backbone radius of gyration as collective variables. For the alanine hexamer, it can be seen from the contour plot, that the method (called d-AFED for "driven" AFED) is able to find the alpha-helical conformation in a simulation of just a few nanoseconds, generating a free energy surface with a deep, well-defined minimum at this conformation.
Recently, we introduced a constant-pressure version of the adiabatic free energy dynamics method in which the cell matrix is used as a set of collective variables as a means of exploring polymorphism in organic molecular crystals. The cell matrix is a useful set of collective variables when the molecules are not polar and interact weakly. This approach, termed Crystal-AFED [see Yu and Tuckerman Phys. Rev. Lett. 107 015701 (2011)], can be combined with the driven AFED scheme discussed above to include internal order parameters as an additional set of collective variables. Using this approach, we are able to predict the free energy differences between the polymorphs of the benzene crystal using simulation lengths of just 5 ns. Shorter simulations (roughly 500 ps) are sufficient to obtain all of the polymorphs. The animation below shows a typical Crystal-AFED trajectory. The figure below the movie shows a plot of the free energies of each polymoprh. Using these free energies, we were able to resolve an ongoing controversy concerning the structure of the benzene II crystal. The figure shows that, although this crystal form has a low lattice energy, it has a high free energy. A benzene II/benzene III mixed structure turns out to have a very low free energy and is likely what is seen in experiments as opposed to benzene II, which is generally not seen.





Software development.

The wide range of projects described above would not be possible without a multi-purpose, flexible, easily extensible software platform in which new developments can be easily implemented and tested and which is also capable of carrying out large-scale simulations on a wide variety of platforms. We have been heavily involved in the development of such a platform (PINY_MD), which is a parallel, object-oriented style software package capable of carrying out force field based, ab initio, and QM/MM calculations (including nuclear quantum effects via path integrals) using state-of-the-art algorithms. The code is freely available under the CPL public license agreement and can be downloaded here. We are currently developing a massively parallel version of the code, called OpenAtom, that employs the charm++ parallelization approach. This code is also now freely available.