Thermostatting in simulation of Electrospray Ionization
An important consideration in MD simulation of ESI droplets is evaporative cooling. Kinetic energy carried away from the droplet by vaporizing solvent molecules lowers the canonical kinetic energy of the droplet bit-by-bit, cooling the droplet. If this continues, the droplet may attain a frozen state where evaporation slows or even stops. Clearly, this evaporative cooling has an effect on the physical properties of the droplet; examples include viscosity and self-diffusion. In an actual ESI apparatus, the heating in the ion source region prevents this evaporative cooling. However, most thermostats employed in MD simulations are unable to capture this effect over the two different phases: droplet (solution) and vapor (expelled solvent and/or solute molecules). See Figure 1. The Lowe-Andersen thermostat is able to keep both phases at a stable temperature, but this thermostat is not implemented in all widely-available MD packages.
Even within MD packages that include it, the Lowe-Andersen thermostat is not currently GPU–accelerated. Thus, if a simulation is to take advantage of GPU-acceleration (obviously desirable for large systems including proteins), simulations must be periodically stopped and velocity distributions must be re-applied which follow the desired temperature. This means that the droplet oscillates between periods of gradual cooling and sudden heating over the course of a full, 1-100ns simulation (Figure 2). Even with a small (3-5%) drop in temperature, it is not clear how these cooling and sudden heating cycles affect the overall evaporation process. It is possible however to observe the evaporation rate using different thermostats. See Figure 3. While a Berendsen thermostat (with or without trajectory stitching) yields a higher evaporation rate than no thermostat at all, a system coupled to the Lowe-Andersen thermostat yields a nearly linear evaporation rate.
The Lowe-Andersen thermostat clearly provides a uniform temperature for the entire system,
but the computational penalty can be high for large systems. Many current ESI evaporation simulations are commonly run with periodic boundary conditions (PBC) for a number of reasons. First and foremost is that PBC simulation code has been GPU-accelerated, and much shorter run times for the same system size is hard to ignore. Care must be taken, however to properly account for and adjust force-field cutoff distances. Particle-mesh Ewald summations must be used to allow for the simulation of non-zero net charged systems. Simulating charged systems without PBC and PME, and without cutoffs (in other words infinite interaction lengths) would be the most realistic choices, but are not always attainable in practice.