Dual-resolution simulations with LAMMPS

So what do I mean by dual-resolution simulations?

One of the most important choices you have to make when trying to model a physical system on a computer is deciding what level of detail is important. Fundamentally, everything may be made up of quarks, gluons and a zoo of fundamental particles (or even vibrations of p-dimensional quantum strings!), but to model the properties of steel girders to calculate the safe loading on a bridge it is good enough (we hope) to treat the steel as a continuous bulk, with properties like compressibility, thermal expansion, and a breaking strain. These are two ends of the heirarchy of multiscale modelling shown in the picture above, but for many biomolecular and materials modelling applications the choice is not so clear-cut.

For example, to simulate chemical reactions of molecules on crystalline surfaces a model which includes charge transfer and chemical bonding such as Density Functional Theory is needed, but this would be far too expensive for all of the 1000s of atoms in the simulation. For these type of applications a combination of DFT and Molecular Mechanics called QM/MM is popular.  Similarly for modelling the dynamics of large biomolecular systems such as proteins embedded in lipid bilayers (the molecular structures that make up animal and plant cell walls), it might be appropriate to model the 'interesting' part of the simulation in atomistic detail, but the 'bulk', or surrounding structure in a more coarse-grained way. These are the dual-resolution models referred to in the title, and in particular one developed by Jon's group called ELBA.

ELBA (ELectrostatics-BAsed) is a model for coarse-grained simulations of lipid systems like the phospholipid above - where groups of around 3-10 atoms are represented as a single spherical bead, and interact with each other via bonded as well as electrostatic (charge and dipolar) forces. Unlike other popular models, it is also easy to introduce atomistic detail where needed, using standard molecular mechanics forcefields like Amber or CHARMM. This is a best-of-both-worlds approach, where the level-of-detail vs performance trade-off can be set by the user.

The main downside to ELBA was that it was only implemented in a relatively unknown code called BRAHMS - ideal for developing the model, but lacking parallelisation and performance optimisations that would allow it to simulate very large systems using HPC. To overcome this, the aim of the eCSE project was to implement ELBA into a mainstream Molecular Dynamics package - in this case LAMMPS, one of the top 10 codes on ARCHER.

The required force calculation code already existed in LAMMPS, so the focus was on implementing a reliable integrator - the code that calculates how the atoms and beads move from one time step to the next. I implemented an algorithm called DLM (named after the authors Dullweber, Leimkuhler and McLachlan) which not only allows for stable long-running calculations with larger timesteps than the in-built LAMMPS integrator, but also supports a range of thermostats and barostats to allow for more realistic simulations than is possible in BRAHMS.

While the dual-resolution approach helps save computational effort it does cause problems for parallel scaling due to load imbalance. The atomistic particles are more expensive to compute so the load balancing code in LAMMPS (which treats all particles as equal) does not work effectively, limiting scalability.  To overcome this I implemented a weighting scheme, where atoms which are more expensive to compute are taken into account by the load balancer. Even for fairly modest simulations - a model of the Bovine Pancreatic Trypsin Inhibitor (BPTI), with 882 atoms surrounded by 6136 ELBA water beads, using 3 ARCHER nodes - I found the weighting scheme to more than double the speed of calculation!

All of this functionality has just been released in the development version of LAMMPS (27Sep16 version), and the DLM integrator is already in the most recent stable release.