I spent the first couple weeks going through tutorials to learn how to use gromacs through tutorials, and now I’m ready to begin my MD simulations on Nafion. For the tutorials I had gone through for GROMACS, I spent most of my time looking at how to perform the actual simulation step and it gave me all of the files beforehand. However, in order to begin simulating Nafion, I needed files specifying the structure and topology of a Nafion polymer chain.
After many discussions with Akhilesh, one of my mentors, I was able to begin building the structure and topology input files for GROMACS. First, I downloaded and installed Avogadro, a program that allows me to build and edit molecules. Using Avogadro, I was able to build a Nafion polymer chain consisting of 10 monomers according to the specified structure below.
Chemical structure of Nafion, https://doi.org/10.1021/jp0761057
A snapshot of 3 of the monomers in the Nafion chain I constructed using Avogadro.
However, this polymer was essentially a linear chain which is unlikely to occur naturally. In an attempt to make things slightly more realistic, I created alternative structures by manually manipulating the polymer into other orientations as shown below.
It was sometimes a little frustrating working with large molecules such as Nafion on Avogadro because Avogadro has a tendency to crash when many atoms exist, but with patience and perseverance I was able to obtain the molecular structures of various nafion chains in the form of ‘.pdb’ files.
With the structure file, I then ran a script, called MKTOP, that uses the GROMACS repositories along with the ‘.pdb’ file to create a topology file that contains parameters for the atoms, bonds, bond angles, and dihedrals in the molecule. All of these parameters are specific to the OPLS-AA (all-atom optimized potentials for liquid simulations) force field being used, however, parameters specific to the Nafion sidechain exist that will allow for a greater accuracy compared to using the OPLS-AA force field.
To implement the new force field parameters, I went through the topology file and replaced the old parameters with the new parameters for only atoms that are a part of the sidechains. I also did this for the bonds, bond angles, and dihedrals of the sidechains. As you can imagine for a molecule with 682 atoms, this took quite some time as the topology file is 6571 lines. In order to be efficient and to save me some sanity, I learned tons of new functionalities for programs like sed and vim to help me edit the topology file without having to go through each line one by one. For example, I can use sed as a way to find and replace certain patterns, but there is also an added functionality that allows me to first search for specific lines and then use the find and replace function only on those lines. Additionally, I learned about and implemented a way in vim to search for specific lines and then fold all other lines to allow me to easily edit only the lines I need. Using a combination of these tools I was able to edit the topology file, and I was able to learn a ton of new text editing tools that I know will make things easier in the future.
When double checking the net charge of the chain with the new parameters, we found that it is not -10. The net charge is expected to be -10 because there are ten sulfonate groups on our nafion polymer that each have a charge of -1. But after further analysis, we found the net charge of the polymer chain to actually be -6.4. In order to get -10, we decided to distribute a charge of -0.1154 to each fluorine on the backbone to achieve this. Because the inconsistency in charge is a result of the new force field, we will run a simulation using purely the OPLS-AA force field to ensure the results are similar and that distributing the remaining charge onto the backbone fluorines will not disrupt the results.
I am looking forward to finally beginning the MD simulations on Nafion in the upcoming weeks!