Earth Sciences & Env.
Engineering & Tech.
Information & Comm. Tech.
Life Sciences & Biotech
Use this search facility to find out more about the profile of our HPC-Europa2 visitors, the type of work they have been doing, and their project achievements.
The main objective of our study was to establish the mechanism of the C-C bond formation step of the self-aldol condensation of propanal reaction catalyzed by Pro-NH2 in water. Such C-C bond formation generates two stereocenters. This means that there are four reactions leading to the four possible stereoisomers; the present project was amied at studying the path yielding the most stable product.
Given the small amount of computational time allocated to the project we could not study all the possible reactions taking place in solution. Nevertheless we have refined an optimal approach to obtain a reliable free energy profile for the C-C bond formation. The blue-moon ensemble thermodynamic integration technique was used in conjuction with Car-Parrinello Molecular Dynamics. It consists in running many simulations at different values of one constraint (in our case C-C distance). The integral of the mean force averaged over the whole simulation gives the free energy for the reaction. With our simulations we were able to study the addition of the propanaldehyde to the enamine to yield the self-aldol condensation product. The barrier is ca. 10 kJ/mol. We have noticed that water plays a fundamental role in stabilizing the transition state for the reaction. This is probably due to stabilization of excess negatve charge via hydrogen bonding. Further studies are needed in order to assess the differences in free energy for the remaining 3 reactions leading to the formation of different stereoisomers.
In this collaborative project we aimed to link the harmonic quantum transition state theory (hQTST) code of Jonsson and Arnaldsson (1) to the DL-FIND code of Kaestner, (2) for use within ChemShell (3). We first aimed to get the rudimentary interface established, before attempting to improve the existing algorithms and parallelising the code. The goal was to yield a versatile programme package that can calculate accurate rates where tunnelling effects are important, for a variety of interesting problems in astrochemistry, biochemistry and surface science.
With the hQTST approximation rates for chemical reactions where quantum mechanical tunnelling through barriers is the dominant mechanism can be calculated quite accurately (1). The DL-FIND code is distributed as a stand-alone optimiser (2), as well as with the ChemShell package (3), which is a flexible, modular code which is interfaced to a variety of quantum chemistry and molecular mechanics (MM) codes and provides various hybrid, cost-efficient schemes applicable to biochemical, material science, and surface science problems.
Currently the hQTST only works with analytical potential energy surfaces, but in principle could work with any routine that yields forces and energies as a function of coordinates. This is an ideal linkage point for DL-FIND/ChemShell, where the energy/force expression can be handled by any quantum mechanics (QM) code, MM-code or with a combined QM/MM approach, where the reactive centre is treated with accurate QM-methods, while the effect of the surroundings is modeled with a computationally less-demanding MM approach. Therefore, the synthesis of hQTST with DL-FIND will yield a versatile tool to study reaction rates in the tunnelling regime. Tunnelling will likely play a role for activated reactions in the interstellar medium where temperatures are very low, but also even at room temperature for various enzymatic and surface reactions where hydrogen atoms are transferred over short distances.
(1) A. Arnaldsson, PhD thesis, Univ. Washington (2007)
We have established the linking of the hQTST code with DL-FIND. We have fully tested this link for using the MNDO code for the energy/force evaluation, and applied this to tunnelling in the malonaldehyde tautomerisation reaction with the cheap, semi-emperical OM1 and AM1 Hamiltonians within MNDO and the B3LYP Hamiltonian with the Gaussian Code.
HQTST aims to find a maximum on the minimum action path from reactant to product. This problem is cast in a Feynman path integral approach, where the classical TS is turned into a discretised closed Feynman path. This translates into finding a first-order saddle point on a potential energy surface of dimension (number of images * degrees of freedom of each image). Originally, this search was performed with a Lanczos algorithm, where the lowest eigenmode is sought for every search step, followed by a conjugate gradient line search along the effect force, where the force along the lowest mode is inverted. It was found that with a reasonable guess for the eigenmode (e.g. the unstable mode of the classical TS, or the unstable mode of a previous converged instanton with fewer images and/or at nearby temperature), the mode does not need to be updated so often. If we update the eigenmode only every 20 steps, 20-70% of computational effort can be saved. A further 50% reduction in computational effort was achieved by using the L-BFGS optimiser instead of conjugate gradients.
It transpired that convergence of the reaction rates at temperatures well below the cross-over temperature (where tunnelling becomes by far the dominant mechanism and the quantum TS delocalises more and more) was slow with respect to the number of images of the discretised Feynman path.
We are currently further investigating how to overcome this, as it is impractical to resort to >100 images for similar systems where one needs a more accurate (and expensive) Hamiltonian. Also, the calculation of the pre-factor of the reaction rates becomes impractical for this large number of images. Since the instanton path itself seems to converge already with fewer images, we are working the implementation of a suitable interpolation scheme for the instanton path as well as calculating the pre-factor from a more modest number of images. First test cases suggest that should be possible, but more extensive testing is necessary.
We have not yet tested the code with a QM/MM example or with different QM-codes, because we did not have any other QM-codes available during this project, but we will soon test this on other compute systems available to us. In the future, we also intend to impement and test the parallelisation of the code by distributing the images among the available nodes. The code can be further parallelised by the intrinsic parallelisation of the QM- and MM-codes that are used.
The objectives of this research sojourn to CINECA were part of a larger molecular dynamics (MD) study of thermotropic carbohydrate and cyclitol liquid crystals. The specific works mainly focussed on the adaption of a force field description for these compounds with regard to molecular dynamics (MD) simulations using GROMACS.
GROMACS, its abilities and force field description format should be explored and understood. Furthermore the process and all necessary components of preparing a simulation for this package should be recognized. Additionally, two different force fields (AMBER and OPLS) and their functional description should be approached with respect to their differences and their applicability to our sugar based compounds.
A MD simulation cannot be run without a thorough preparation which includes numerous steps. The force field parameter containing topology has to be set up. A reasonable ensemble of atom positions, the initial configuration, has to be generated. Finally, the parameters for the simulation run have to be carefully chosen. The experience of Prof. Zannoni’s group should ease this process and avoid starting from scratch.
The force fields AMBER and OPLS were evaluated for our purpose, which are related in their functional form. Both have a similar way of describing bonded and non-bonded interactions. Bonds and angles are represented by harmonic potentials, dihedrals by Fourier series or cosine sums, respectively, and the nonbonded interactions (Coulomb and van der Waals) via a Lennard-Jones potential.
A major difference between them is the assignment of the partial charges. While the OPLS charges are applied in a generalized fashion, while the AMBER force field approaches the partial charges based on a quantum mechanical population analysis.
Preparing a MD simulation in GROMACS consists of three aspects, the topology, the configuration, and the simulation parameters.
The force field description of a chemical compound is created, in GROMACS terms, as a topology file. It contains the definitions of atom-, bond-, angle-, and dihedral types. The atom type consists of the default partial charges, e and s for the Lennard-Jones potential, and a bond type name. The bond-, angle-, and dihedral types are defined with these bond type names of the participating atoms and an integer code for the functional description used in GROMACS. Sometimes it is possible to derive force field parameters from literature; however, their dimension units have to be considered carefully.
Along with the type definitions, the actual atoms, bonds, angles, dihedrals, and pairs for 1-4 interactions (intramolecular, non-bonded interactions between atoms, that are separated by two connecting atoms) are specified. Each atom in a molecule, identified by a unique atom name, is assigned an atom type. The partial charges in a molecule are grouped into neutral charge groups. Bonds, angles, dihedrals, and 1-4 pairs have to be stated explicitly. The relation to the type descriptions mentioned above is implicitly created by the bond type name of the assigned atom types.
An important aspect of the simulation of order evolution is to avoid pre-ordered starting configurations by stacking tools. The PACKMOL tool enables the generation of totally disordered configurations. With these configurations stepwise cooling simulations can be performed to study the onset of order.
Due to the manifold options that are possible in MD simulations a reasonable choice was made with the experience of the Prof. Zannoni’s group. To keep the temperature constant the velocity rescale algorithm was chosen, that adds a statistical term to the Berendsen thermostat. Isotropic pressure coupling was obtained by the Berendsen barostat. The Particle Mesh Ewald (PME) method was applied to evaluate electrostatic long range influences. A cut-off algorithm mimics the van der Waals short range interactions.
The built topology, configuration, and simulation parameter files are compiled with the GROMACS precompiler “grompp” into a final simulation input file.
The actual MD application “mdrun” performs the simulation on this file. The workflow consists of the preparation of a starting configuration, a minimization using the same constraints, neighbour searching, electrostatics, and VDW parameters as in the actual production run. A short equilibration run of a few nanoseconds is performed until observables like potential energy, volume, and density remain stable.
The production run with progressive cooling uses a stepwise cooling schedule to sample every temperature for a sufficient time.
GROMACS offers a wide variety of analysis tools. They operate mainly on the energy or trajectory output files of “mdrun”. Trajectories are positions stored over time. The new tool developed by my co-worker Mr Breuers could be used to analyse the trajectories with respect to their orientational order parameter.