[Blog #1] My MD simulations are too hot to handle, at least for my own computer

I have now finished the first few weeks of my research internship at LBNL! So far, I have been on a non-stop learning journey, building up my background knowledge and my skills to reach the point where I can efficiently perform the tasks needed for the research project.


As much of my research involves studying the properties of Nafion, a PFSA ionomer used as the proton exchange membrane in fuel cells, Dr. Su sent me a collection of research papers to help me gain a deeper understanding of some concepts central to the project, including the morphology of PFSA ionomers, the effects of different ions in ionomers, and the use of x-ray scattering and x-ray absorption spectroscopy in studying ionomers. Originally, my research project was to study the effects of different counterions in Nafion using x-ray scattering methods, however I am unable to perform the necessary experiments in the lab as the internship is remote. As a way around these limitations, we decided to use molecular dynamics (MD) methodologies to simulate Nafion and study the effects of different counterions all on a computer. One of the main hopes we have is to eventually extract the necessary spectroscopy data from the MD simulations and compare them to some of the experimental data already obtained.


Monomer of Nafion, 10.4028/www.scientific.net/JNanoR.5.31

In order to get started with the MD simulations, Dr. Su and I met with Akhilesh Paspureddi, a Ph.D. student at UT Austin who has experience simulating Nafion using molecular dynamics. Since Akhilesh is very experienced with the GROMACS software package for MD simulations, I have spent the first few weeks becoming comfortable with GROMACS learning from both tutorials, as well as from Akhilesh himself.


GROMACS was developed by the Department of Biophysical Chemistry at the University of Groningen, and it essentially functions by simulating the Newtonian equations of motion for very large and complex systems containing up to hundreds of millions of molecules. I learned that three files are required for GROMACS: a structure file, a topology file, and a parameters file. The structure file specifies the coordinates of atomic sites. The topology file specifies the forcefield, which holds information on both the bonded and nonbonded atomic interactions. Lastly, the parameters file holds details for the simulation itself, such as how long the simulation will be run and what algorithms will be used. In order to perform a simulation itself, there are three major steps. The first step is a minimization step to minimize the potential energy in the system and ensure the molecules are in an appropriate starting position. The second step is an equilibration step to stabilize both the temperature and the pressure of the system. The third step is the production step where the actual simulation is run.


To get started with GROMACS, I first completed the tutorials developed by the CNPS Lab at the University of Rijeka that walked through simulating simple molecules like water and methane.


 I created a methane molecule solvated in a 3.2 x 3.2 x 3.2 nm box of water.


I then minimized the potential energy of the system. The graph below is a plot showing the potential energy converging as the number of steps of the minimization algorithm increased.



I proceeded to perform an equilibration to stabilize the temperature of the system as the left plot below of temperature vs. time shows the temperature converging to around 298 K. I then continued with a pressure equilibration. It is a bit difficult to see in the right plot below of pressure vs. time, but the pressure fluctuations are centered closely to 1 bar, which matches expectations.



After completing the minimization, equilibration, and production steps, I was then able analyze specific properties of the system such as the density and the radial distribution function. To do so, GROMACS has built-in functions that allow for a variety of properties to be studied. The calculated density was around 993 kg/m3, as shown by the GROMACS output line below, which is very close to the expected value of 995 kg/m3 for the specific force field used.



Additionally, the radial distribution function was calculated and plotted below, which shows the density of water molecules around the methane atom.



After completing more of these simpler tutorials, I moved to a much larger protein molecule by following the tutorials made by Justin Lemkul, a researcher at Virginia Tech. However, before I started this tutorial, I wanted to avoid using my local computer. Up until this point, I had only run GROMACS using my personal computer, and I found the simulations were becoming pretty expensive. For example, the calculations for ten methane molecules in a box of water took almost a day and a half using 98% of my CPU throughout its duration! With my computer fan being on full blast for days on end combined with the lack of CPU power to use my computer for anything else, I decided it was time to figure out how to use the Lawrencium cluster we have available at LBNL to do these calculations. A cluster is essentially a bunch of computers combined in a way where they can interact with each other and function in parallel to one another. So instead of being limited to the 4 processing cores I have on my personal computer, I can use a ton more cores on the cluster (I typically use 72), which reduces the calculation time from over a day down to a couple hours!

I originally thought this would not be too difficult to get started on because of my experience using the cluster for DFT calculations, but it did take a bit of time to figure out a working job script for using GROMACS on the cluster. After a couple hours of reading more about mpi processing, Slurm Workload Manager, and details about using GROMACS on a cluster, I was able to gain a better understanding of how to submit jobs, and I was finally able to figure out a working job script through trial and error. Now I could simulate the bigger molecules, such as the lysozyme below, using the same methodology used to simulate methane in water!


Visualization of the hen egg white lysozyme solvated in a box of water. The lysozyme is the blue patch in the middle surrounded by the red water molecules.


As you can see, there is substantially more molecules here, so I am very thankful I was able to use the cluster for this tutorial. As I am now much more comfortable with GROMACS, I am ready to work on simulating Nafion in the upcoming weeks!