(This tutorial is related to work published at https://doi.org/10.3390/nano12193379).
First, we will build a cellulose nanofibril structure using the program Cellulose-builder. The program offers a user-friendly interface and a set of bash scripts to build the cellulose structure in PDB format with different cellobiose chain lengths and geometries (crystalline, fibril, and monolayer). The program allows imposing periodicity along the chain axis for all model types and 3-dimensional periodicity for the crystalline model. The program supports geometry building for various cellulose polymorphs (I-α, I-β, II, and III(I)).
We will build a cellulose nanofibril with 10 cellobiose units of the 1-4 linked I-β polymorph for this tutorial.
Cellulose builder web interface:
If you prefer to avoid the fuss of downloading and running a script, click here to access the online GUI for geometry construction (Go to the next subsection for a command-line bash script version for structure building). Below is an image of the webserver with queries chosen to build a single nanofibril of 1-4 linked I-β cellulose with 10 cellobiose units (20 D-glucose molecules) per chain. We want to model to represent a slice of the cellulose fibre, so we choose to link the terminal glucose molecules.
Select Input: fibril
Input 1: 10
NOTE: When choosing to build a crystalline representation of cellulose, it is important to go through the article and check the significance of the three input integer options, which define the number of chains on origin and centre layers, and the chain length. Below is an example of parameters to build a structure with 6 origin layers, 4 chains/origin centre layers, and 5 cellobiose units/chains.
After clicking on the submit button, you can download the results in a zip file that includes the structure output in PDB format (crystal.pdb) and other parameter files supporting the AMBER MD package (not significant for this tutorial). Navigate to the download directory using the terminal and use the following command to extract results from the compressed file:
Cellulose Builder command-line scripts:
If you prefer to dive into the scripts and look at how the geometries and different structure types are built, this is the step for you. The bash script source code of cellulose builder can be downloaded here (GNU license). Using the scripts to build the model has few installation dependencies (Octave, VMD and psfgen). Check the official documentation of the program here for instructions on how to install the pre-requisites. Remember to add the dependencies executables to the PATH to your “bashrc” or “bash-profile”. Okay, now that we have all the dependencies installed, let’s extract the scripts from the compressed file using the command-line terminal.
For tarball file: tar -xzf cellulose-builder_July_2013.tar.gz
For zip file: unzip cellulose-builder_July_2013.zip
Next, start by navigating to the cellulose-builder_July_2013 directory. We first modify the file input.inp to choose the I-beta polymorph and chain periodicity. (You can also use any text editor of your choice)
vi input.inp (below image shows the choice of parameters for I-beta fibril with periodicity along the axis of chain)
Now run the bash script using the following command:
./cellulose-builder fibril 10
If all dependencies and executable paths are correctly defined, you will see a new directory named crystal, inside which there will be a PDB format structure file named crystal.pdb.
Now let’s take a closer look at the structure to get an idea of the orientation and interactions.
You can navigate through the structure to look at the orientation and intramolecular interaction within the chain and across parallel chains.
Now let’s set up the first MD simulations for the fibril system. A trivial step here is to generate the Gromacs compatible topology files for the cellulose system. Here we use a single D-glucose molecule parameterization to describe the bonded and non-bonded parameters of the system. The parameters (bonds, angles, dihedrals and partial charges) have been already derived with methods compatible and validated for Amber14 force fields, so we will not go into those detail in this tutorial. We will use the Amber14 force field parameters for this tutorial, which can be downloaded as amber14-mfc-pp.tgz. Unarchive the force field tarball into a new directory:
mv ~/Downloads/amber14-mfc-pp.tgz force-fields/
tar -xzf amber14-mfc-pp.tgz
There should now be a subdirectory named amber14sb_OL15.ff and two files, residuetypes.dat and specbond.dat.
NOTE: If you have time, have a look at the files inside the amber14sb_OL15.ff directory. The glucose molecule parameters are described in the following files:
- New-glc.rtp: Residue topology with partial charges and bond descriptions
- New-glc.r2b: Residue name translation database
- New-glc.arn: Atom renaming database
For a detailed description of the file types, check this link on the Gromacs webpage.
We want the Gromacs library to access these force fields from any working directory of the computer, so let’s add the path to the bashrc file of your system.
echo export GMXLIB=”<path-to-force-field-directory>” >> ~/.bashrc
or run this command to add the path to the environment variables.
Now, let’s start generating the topologies for a cellulose nanofibril in water. Please create a new directory, cellullose-fibril-water and copy the structure file we generate (crystal.pdb) to this new directory. Under this directory, we set up our MD simulations using an approach that allows linking the terminal residues of the cellulose chain across the periodic box.
Prepare the PDB file to match Gromacs RTP description
sed -i "s/BGLC/IBG /" <structure-filename>.pdb
Generate preliminary topology files to be used for the [atom] section description
gmx_mpi pdb2gmx -f crystal.pdb -o pdb2gmx.gro -p topol-pdb2gmx.top -water tip3p
The pdb2gmx tool will process the structure and prompt to choose a force field:
(Choose option 1:AMBER14sb-OL15 + beta-glucose GLYCAM06h)
Ignore the long bond warning here since pdb2gmx does not support periodic covalent bonding, and we will re-define these using x2top a few steps later.
Important: Compute the cellulose fibril chain length from the monomer count to define the length of box dimension "c" to allow covalent bonding of the terminal residues. In our case, we have 10 monomer/chain ==> 10.38 nm length.
gmx_mpi editconf -f pdb2gmx.gro -o mfc-20mono.gro -box 7 8 10.38
Now let's generate a topology with periodic bonded parameters using "x2top" instead of pdb2gmx.
gmx x2top -f mfc-20mono.gro -o topol-x2top.top -ff select -pbc -name IBG -noparam -alldih
(Choose option 1:AMBER14sb-OL15 + beta-glucose GLYCAM06h)
The next step is to add some corrections to the topology files and generate the final cellulose topology file "newtopol.top".
sed -i '/\[ dihedrals \]/,/\[ system \]/s/1 *$/9/' topol-x2top.top
sed -n '/Include forcefield parameters/,/\[ bonds \]/p' topol-pdb2gmx.top | sed '/bond/d' >& newtopol.top && sed -n '/\[ bonds \]/,/\[ system \]/p' topol-x2top.top | sed '/\[ system \]/d' >> newtopol.top && sed -n '/Include Position restraint file/,//p' topol-pdb2gmx.top >> newtopol.top
The rest of the MD workflow is now similar to the other Gromacs tutorials. We fill the fibril cell with water:
gmx_mpi solvate -cp mfc-20mono.gro -cs spc216.gro -p newtopol.top -o box-sol.gro
The solvated system is now assembled, and the output structure file box-sol.gro can be visualized with PyMol. Also, check for the modifications in the [ molecules ] section of the topology file (newtopol.top), where the number of waters should be included.
Before we begin the production dynamics, we want to relax the system and remove any steric clashes, using an energy minimization run. A key parameter necessary to simulate periodic/infinite cellulose chain parameters is the addition of "periodic_molecules" to the mdp file. Generate the binary input tpr file using mdp parameter file em.mdp
gmx_mpi grompp -f em.mdp -p newtopol.top -c box-sol.gro -o em
We can now proceed with the relaxation MD step. Try to run the next step on the cluster, as it may take about an hour on a regular computer. You will have to copy all the directories we created to the cluster (force fields and cellulose-fibril-water). Important, do not forget to add the line export GMXLIB=”<path-to-force-field-directory>” to your bashrc file.
gmx_mpi mdrun -deffnm em
Evaluate if your minimization run was successful by checking the em.log file. If not, then try to check what is happening based on the log output. Consider doing a quick analysis of the evolution of the potential energy of the system.
gmx_mpi energy -f em.edr -o em.xvg
Choose option number for Potential and then “0” to terminate the analysis input. Use Xmgrace or any other program of choice to visualize the energy plot.
You can also visualize the optimized structure in the em.gro file. But before that, we will fix any periodicity issues in the structure since Gromacs may not always put back the structure as a whole after implementing a periodic box.
echo “Cellulose System” | gmx_mpi trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro
A successful run will ensure we have a good starting model for the equilibration step. We will perform a two-step equilibration of the solvated fibril, starting with temperature equilibration with the canonical (NVT) ensemble and followed by an isothermal-isobaric ensemble (NPT) of 500ps each. The short time is sufficient for equilibrating small/medium systems, but we check the output to re-affirm the plan.
We now repeat the grompp and mdrun commands using a different mdp file nvt.mdp:
gmx_mpi grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt
gmx_mpi mdrun -deffnm nvt -dlb yes
At the end of a successful run, analyze the temperature evolution and output structure:
gmx_mpi energy -f nvt.edr -o nvt.xvg
Choose the option "Temperature" followed by "0" to end the input prompt. Now, using any plotting tool of choice (Xmgrace), you can visualize the temperature fluctuation.
Echo “Cellulose System” | gmx_mpi trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro
After bringing the system's temperature to a stable range of around 300K, we will now use a barostat in combination with the thermostat to model the system with an environment that closely resembles the experimental conditions. In the second step equilibration, we will run the NPT ensemble run for 500 ps with a V-rescale thermostat, weakly-coupled Berendsen barostat and position restraints for heavy atoms. Note that we implement semi-isotropic pressure coupling to separate the way the fibril feels pressure along the chain axis (Z) and the other two axes. The pressure coupling is isotropic in the X and Y direction. The NPT run parameter file is available as npt.mdp.
; Pressure coupling
pcoupl = Berendsen ;Parrinello-Rahman for production run (NPT)
pcoupltype = semiisotropic ; uniform scaling of box vectors in x-y and independent z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
Now, let’s start the pressure equilibration run.
gmx_mpi grompp -f npt.mdp -p newtopol.top -c nvt.gro -r nvt.gro -t nvt.cpt -o npt
gmx_mpi mdrun -deffnm npt -dlb yes
Analyze the pressure equilibration using the following command:
gmx_mpi energy -f npt.edr -o npt.xvg
Choose the options "Pressure" followed by "0" to end the input prompt and visualize the plot in the output file npt.xvg
Next, we again make the structure whole using the trjconv tool:
echo “Cellulose System” | gmx_mpi trjconv -f npt.gro -s npt.tpr -pbc mol -center -o npt.gro
As expected, the pressure value will fluctuate widely, but the average should be closer to the reference pressure value of 1 bar. The reason behind the fluctuations has been explained in detail in the Gromacs manual and tutorial.
We use the same NPT ensemble for the production MD run but without any position restraints and a tighter pressure coupling for using the Parrionello-Rahman barostat. You can either copy the npt.mdp to md.mdp and make the below changes or download the parameter file md.mdp.
pcoupl = Parrinello-Rahman
Production run & Analysis
Okay, finally time to start the production run.
gmx_mpi grompp -f md.mdp -p newtopol.top -c npt.gro -o topol
gmx_mpi mdrun -s topol -dlb yes
Make sure to check the run ends successfully. In case of time limit issues on the cluster, it is possible to restart the MD run from the last saved time-point (state.cpt) using the following command:
gmx_mpi mdrun -s topol -dlb yes -cpi state.cpt
Once completed, make the output structure in the file confout.gro whole again using trjconv tool. Check the structure with PyMol for any visually distinct feature in the fibril model.
echo “Cellulose System” | gmx_mpi trjconv -f confout.gro -s topol.tpr -pbc mol -center -o confout-10ns.gro
In this section, we only analyse the stability of the model during the 10 ns run. It is always a good idea to evaluate the RMSD fluctuation of the most relevant part of the system (cellulose fibril in our case).
echo “Cellulose Cellulose” | gmx_mpi rms -f traj.trr -s topol.tpr -o rmsd
In the plot, there may be an initial rise in root-mean-square fluctuations but it should stabilize within 2 ns with only small fluctuations. Next, we will check the potential energy evolution over the trajectory if the stability in the RMSD curve matches the energy landscape.
gmx_mpi energy -f ener.edr -o energy.xvg
Choose the options "Potential", "Temperature", and "Density" followed by "0" to end the prompt. Now open the plot file and evaluate if the trajectory needs to be extended based on the energy landscape and density fluctuations.
Click here to move to the tutorial for building a polypropylene solvent box