Steering in computational science: mesoscale modelling and simulation

This paper outlines the benefits of computational steering for high performance computing applications. Lattice-Boltzmann mesoscale fluid simulations of binary and ternary amphiphilic fluids in two and three dimensions are used to illustrate the substantial improvements which computational steering offers in terms of resource efficiency and time to discover new physics. We discuss details of our current steering implementations and describe their future outlook with the advent of computational grids.Comment: 40 pages, 11 figures. Accepted for publication in Contemporary Physic

Computer simulations of polymers and gels

Computer simulations have become a vital tool in modern science. The ability to reliably move beyond the capabilities of experiment has allowed great insights into the nature of matter. To enable the study of a wide range of systems and properties a plethora of simulation techniques have been developed and refined, allowing many aspects of complex systems to be demystified. I have used a range of these to study a variety of systems, utilising the latest technology in high performance computing (HPC) and novel, nanoscale models. Monte Carlo (MC) simulation is a commonly used method to study the properties of system using statistical mechanics and I have made use of it in published work [1] to study the properties of ferrogels in homogeneous magnetic fields using a simple microscopic model. The main phenomena of interest concern the anisotropy and enhancement of the elastic moduli that result from applying uniform magnetic fields before and after the magnetic grains are locked in to the polymer-gel matrix by cross-linking reactions. The positional organization of the magnetic grains is influenced by the application of a magnetic field during gel formation, leading to a pronounced anisotropy in the mechanical response of the ferrogel to an applied magnetic field. In particular, the elastic moduli can be enhanced to different degrees depending on the mutual orientation of the fields during and after ferrogel formation. Previously, no microscopic models have been produced to shed light on this effect and the main purpose of the work presented here is to illuminate the microscopic behaviour. The model represents ferrogels by ensembles of dipolar spheres dispersed in elastic matrices. Experimental trends are shown to be reflected accurately in the simulations of the microscopic model while shedding light on the microscopic mechanism causing these effects. These mechanisms are shown to be related to the behaviour of the dipoles during the production of the gels and caused by the chaining of dipoles in magnetic fields. Finally, simple relationships between the elastic moduli and the magnetization are proposed. If supplemented by the magnetization curve, these relationships yield the dependencies of the elastic moduli on the applied magnetic field, which are often measured directly in experiments. While MC simulations are useful for statistical studies, it can be difficult to use them to gather information about the dynamics of a system. In this case, Molecular Dynamics (MD) is more widely used. MD generally utilises the classical equations of motion to simulate the evolution of a system. For large systems, which are often of interest, and multi-species polymers, the required computer power still poses a challenge and requires the use of HPC techniques. The most recent development in HPC is the use of Graphical Processing Units (GPU) for the fast solution of data parallel problems. In further published work [2], I have used a bespoke MD code utilising GPU acceleration in order to simulate large systems of block copolymers(BC) in solvent over long timescales. I have studied thin films of BC solutions drying on a flat, smooth surface which requires long timescales due to the ’slow’ nature of the process. BC’s display interesting self-organisation behaviour in bulk solution and near surfaces and have a wide range of potential applications from semi-conductors to self-constructing fabrics. Previous studies have shown some unusual behaviour of PI-PEO diblock co-polymers adsorbing to a freshly cleaved mica surface. These AFM studies showed polymers increasing in height over time and proposed the change of affinity of mica to water and the loss of water layers on the surface as a driver for this change. The MD simulation aimed to illuminate the process involved in this phenomena. The process of evaporation of water layers from a surface was successfully simulated and gave a good indication that the process of solvent evaporation from the surface and the ingress of solvent beneath the adsorbed polymer caused the increase in height seen in experiment

Large-scale parallelised boundary element method electrostatics for biomolecular simulation

Large-scale biomolecular simulations require a model of particle interactions capable of incorporating the behaviour of large numbers of particles over relatively long timescales. If water is modelled as a continuous medium then the most important intermolecular forces between biomolecules can be modelled as long-range electrostatics governed by the Poisson- Boltzmann Equation (PBE). We present a linearised PBE solver called the "Boundary Element Electrostatics Program"(BEEP). BEEP is based on the Boundary Element Method (BEM), in combination with a recently developed O(N) Fast Multipole Method (FMM) algorithm which approximates the far-�field integrals within the BEM, yielding a method which scales linearly with the number of particles. BEEP improves on existing methods by parallelising the underlying algorithms for use on modern cluster architectures, as well as taking advantage of recent progress in the �field of GPGPU (General Purpose GPU) Programming, to exploit the highly parallel nature of graphics cards. We found the stability and numerical accuracy of the BEM/FMM method to be highly dependent on the choice of surface representation and integration method. For real proteins we demonstrate the critical level of surface detail required to produce converged electrostatic solvation energies, and introduce a curved surface representation based on Point-Normal G1-continuous triangles which we �find generally improves numerical stability compared to a simpler surface constructed from planar triangles. Despite our improvements upon existing BEM methods, we �find that it is not possible to directly integrate BEM surface solutions to obtain intermolecular electrostatic forces. It is, however, practicable to use the total electrostatic solvation energy calculated by BEEP to drive a Monte-Carlo simulation

Simulating molecular docking with haptics

Intermolecular binding underlies various metabolic and regulatory processes of the cell, and the therapeutic and pharmacological properties of drugs. Molecular docking systems model and simulate these interactions in silico and allow the study of the binding process. In molecular docking, haptics enables the user to sense the interaction forces and intervene cognitively in the docking process. Haptics-assisted docking systems provide an immersive virtual docking environment where the user can interact with the molecules, feel the interaction forces using their sense of touch, identify visually the binding site, and guide the molecules to their binding pose. Despite a forty-year research e�ort however, the docking community has been slow to adopt this technology. Proprietary, unreleased software, expensive haptic hardware and limits on processing power are the main reasons for this. Another signi�cant factor is the size of the molecules simulated, limited to small molecules. The focus of the research described in this thesis is the development of an interactive haptics-assisted docking application that addresses the above issues, and enables the rigid docking of very large biomolecules and the study of the underlying interactions. Novel methods for computing the interaction forces of binding on the CPU and GPU, in real-time, have been developed. The force calculation methods proposed here overcome several computational limitations of previous approaches, such as precomputed force grids, and could potentially be used to model molecular exibility at haptic refresh rates. Methods for force scaling, multipoint collision response, and haptic navigation are also reported that address newfound issues, particular to the interactive docking of large systems, e.g. force stability at molecular collision. The i ii result is a haptics-assisted docking application, Haptimol RD, that runs on relatively inexpensive consumer level hardware, (i.e. there is no need for specialized/proprietary hardware)

Ligand recognition by the major urinary protein

Molecular Dynamics (MD) and Quartz Crystal Microbalance (QCM) techniques can provide unique insights into what drives protein-ligand association. The major urinary protein (MUP) binds small ligands in a deeply buried hydrophobic pocket. Detailed calorimetric studies have shown that ligand binding is driven by enthalpic effects, not entropic effects [1]. Previous studies have shown that this is due to 'dewetting' of the binding site cavity even in the absence of ligands, and have also characterised the complex changes in molecular flexibility that accompany ligand binding-features that may be correlated with NMR data [2]. Recent MD revealed the hydration effects of apo-MUP and also shown where certain regions of MUP become more flexible upon ligand binding. They have also shown a water molecule remains close to the tyrosine in the binding pocket [2]. In our current MD studies and OCM experiments we have used wild type and 2 different mutants of MUP to study the binding effects of the ligand IBM. The first mutant has an OH group removed from the binding site of MUP (i.e. tyrosine to phenylalanine (Y120F)). The second mutant has an extra OH group in the binding site (i.e. alanine to serine (A103S)). For all three systems the hydration and flexibility upon ligand binding has been analysed. The hydration analysis from MD reveal (from radial distribution curves and hydration density maps) there is a small density of water that remains even without the presence of the ligand for the WT MUP whereas a larger density of water remains in the binding cavity of the A103S hydrophilic MUP simulation. The results are based on the average structure generated from the 1 mus simulations. The Y120F MUP simulations reveal that there is no water molecules present in the binding cavity. However, as protein molecules are very dynamic in nature, water molecules are observed to hop in and out of the binding pockets for both mutant MUP (but not WT MUP) simulations over the 1 mus simulations. On the other hand the experimental QCM results reveal that on ligand binding no water loss is observed for Y120F mutant MUP whereas A103S and WT MUP have about 2 water molecules which are lost in the binding cavity. The flexibility results from the MD simulations reveal that WT MUP have some residues which increase in flexibility whilst other residues which decrease in flexibility on ligand binding. However, the Y120F hydrophobic MUP show an overall decrease in flexibility whereas the A103S MUP shows an overall increase in flexibility on ligand binding. In contrast the experimental OCM and AFM results reveal that there is an increase in flexibility on ligand binding to all 3 different types of MUP molecules. The experimental and the simulation data have shown a variation in results but it is to be noted that the results cannot be directly compared as the analytical experiments are a surface based techniques whereas the MD simulations do not involve a surface. However, the contrast observed between computer simulation and experiments has revealed important information on the ligand binding effects on MUP. [1] Bingham, R.J., J.B.C. Findlay, S.Y. Hsieh, A.P. Kalverda, A. Kjeliberg, C. Perazzolo, S.E.V. Phillips, K. Seshadri, C.H. Trinh, W. B. TurnbulI, G. Bodenhausen, and S.W. Homans. 2004. Thermodynamics of binding of 2-methoxy-3-lsopropylpyrazlne and 2- methoxy-3-lsobutylpyrazine to the major urinary protein. J. Am. Chem. Soc. 126:1675-1681. [2] Barratt, E., R.J. Bingham. D.J. Warner, C.A. Laughton, S.E.V. Phillips, and S.W. Homans. 2005. Van der Waals interactions dominate ligand-protein association in a protein binding site occluded from solvent water. J. Am. Chem. Soc. 127:11827-11834

Development of High Performance Molecular Dynamics with Application to Multimillion-Atom Biomass Simulations

An understanding of the recalcitrance of plant biomass is important for efficient economic production of biofuel. Lignins are hydrophobic, branched polymers and form a residual barrier to effective hydrolysis of lignocellulosic biomass. Understanding lignin\u27s structure, dynamics and its interaction and binding to cellulose will help with finding more efficient ways to reduce its contribution to the recalcitrance. Molecular dynamics (MD) using the GROMACS software is employed to study these properties in atomic detail. Studying complex, realistic models of pretreated plant cell walls, requires simulations significantly larger than was possible before. The most challenging part of such large simulations is the computation of the electrostatic interaction. As a solution, the reaction-field (RF) method has been shown to give accurate results for lignocellulose systems, as well as good computational efficiency on leadership class supercomputers. The particle-mesh Ewald method has been improved by implementing 2D decomposition and thread level parallelization for molecules not accurately modeled by RF. Other scaling limiting computational components, such as the load balancing and memory requirements, were identified and addressed to allow such large scale simulations for the first time. This work was done with the help of modern software engineering principles, including code-review, continuous integration, and integrated development environments. These methods were adapted to the special requirements for scientific codes. Multiple simulations of lignocellulose were performed. The simulation presented primarily, explains the temperature-dependent structure and dynamics of individual softwood lignin polymers in aqueous solution. With decreasing temperature, the lignins are found to transition from mobile, extended to glassy, compact states. The low-temperature collapse is thermodynamically driven by the increase of the translational entropy and density fluctuations of water molecules removed from the hydration shell