Simulation with induced pressure difference

Now you will come back to the simulation that you submitted at the beginning. First, let's learn something about the background and motivation for this type of simulation.

In experiments, a typical method to study water channels is to set up different concentrations of (impermeable) solutes on the two sides of the channel, thus giving rise to an osmotic pressure difference, and to measure the net directional water flux. Obviously, the relationship between the water flux and the osmotic pressure difference cannot be studied by equilibrium simulations.

It is known that an osmotic pressure is equivalent to a hydrostatic pressure. Therefore, if one can generate a hydrostatic pressure difference in MD simulations, one could mimic the experiments mentioned above. This can be done by applying a constant force along the $z$ direction on a layer of bulk water molecules. The force will induce a pressure gradient in the bulk water layer, which results in a pressure difference on the two sides of the nanotubes due to the periodic boundary conditions (refer to Fig. 2). One appealing feature of this method is that the amount of pressure difference can be easily calculated and tuned, thus the results can be compared to experiments quantitatively.

In the submitted simulation, a constant force of 0.4 Kcal/mol/Å along the $+z$ direction was applied to a 5.4 Å-thick water layer. In every time step of the simulation, the force will be applied to the O atoms of water molecules whose coordinates (in the unit cell of the periodic system) satisfy $z>12.5$Å or $z<-12.5$Å. This is a fairly large force, but it is needed to induce a very fast water flux that can be observed in a very short simulation time (practically you can only afford a 40 ps simulation for the purpose of this tutorial). You may take a few moment to inspect the NAMD configuration file sim_short.conf for this simulation.

1. Check whether the simulation has already finished. If not, wait a few minutes for it to finish (after 40,000 steps).

2. In VMD, create a new molecule with the files cnt.psf and sim_short.dcd (the DCD file generated by your simulation).

3. Play the trajectory and observe the water flow. Since water molecules in the nanotubes move concertedly, you can label any of them in a nanotube, and its movement will indicate the movement of the whole water chain. Did you observe the drift of the water chains toward $+z$ direction in the trajectory? What is the net water flow through the nanotubes during this simulation?

NOTE: The net water flow is in general different from the number of permeation events. For example, when each water molecule moves one step in the single file, effectively a water molecule is transferred, thus resulting in a net water flow of one water molecule, although this move may not cause a permeation event. In the equilibrium simulation you looked at earlier, there was little net water flow, although you observed a lot of permeation events. In contrast, here in the present (short) trajectory, you might not see any individual water molecule crossing all the way through the nanotubes, but you can observe a directional water flow.

The trajectory of your simulation (40 ps) is too short for good statistics. Therefore, here we provide to you another trajectory, which was run at exactly the same conditions as the simulation you ran, except that this one is longer (1 ns), which allows you to observe more water flow.

4. In VMD, delete all the frames of the current molecule (which should be cnt.psf), by going to Molecule $\to$ Delete Frames $\to$ Delete. Then load into the DCD file sim.dcd, by first selecting the cnt.psf molecule and then calling the File $\to$ Load Data Into Molecule... menu item.

5. Observe the water flow in the long trajectory. Since the total water flow is fairly large, it will be too tedious to count it by hand. Therefore, we have provided to you a script, flow.tcl, to count the water flow. In VMD, make sure that the 1 ns trajectory is the ``top'' molecule, and then type in the TKCon window source flow.tcl. The script will report the amount of net water flow and its direction.

NOTE: While the script is running, VMD will freeze and will not respond to other commands. Depending on your machine, this may take 10-20 seconds.

Exercise: Can you calculate the hydrostatic pressure difference in this simulation? The pressure difference is $\Delta P=nf/A$. Here you know that $f$ is 0.4 Kcal/mol/Å, the unit cell of the periodic system has dimensions of 23 Å$\times$ 19.9 Å$\times$ 30.4 Å, and $n$, the number of water molecules in the 5.4 Å-thick layer (water molecules on which you are applying force), can be roughly estimated by the volume of the layer and the molar volume of water (55.5 mol/l). What is your calculated pressure difference? How does it compare to one atmosphere (10$^5$ Pa)? Now you may realize how ``hard'' it is to induce a significant water flow within 40 ps.

6. You may also look at water orientation in the nanotubes as water is being conducted. Does the fast flow alter the water orientation that you had observed in the equilibrium case?