(This tutorial is related to work published at https://doi.org/10.3390/nano12193379).
This tutorial will build a composite model with a hydrophilic cellulose nanofibril enclosed by a hydrophobic polypropylene matrix. In the previous tutorials, we built the structure models for the individual components (cellulose nanofibril and polypropylene matrix). Here we will build the composite model and use the simulated annealing method to generate a computational model with the target densities.
It is recommended to first go through the tutorials of building models for cellulose nanofibril in a water box and a polypropylene box. If you have already completed those or prefer to skip that part then download the geometries for cellulose nanofibril with 10 cellobiose units and for the (poly)propylene chain composed of 10 monomers.
Let's begin with a new directory for the biocomposite system.
mkdir biocomp-mfc20-pp10 && cd biocomp-mfc20-pp10
Use the below commands to copy the nanofibril structure and topologies from the previous tutorial:
cp ~/cellulose-fibril-water/newtopol.top newtopol-mfc20.top && cp ~/cellulose-fibril-water/posre.itp posre-mfc20.itp
sed -i "s/posre/posre-mfc20/ ; /SOL/d" newtopol-mfc20.top
cp ~/cellulose-fibril-water/mfc-20mono.gro mfc-20mono.gro
Next, copy the polypropylene single-chain structure and topologies to the current directory
cp ~/PP-10mono/single-chain-10mono.gro pp-single-chain-10mono.gro
cp ~/PP-10mono/single-chain-10mono.itp pp-single-chain-10mono.itp
cp ~/PP-10mono/single-chain-10mono.gro pp-single-chain-10mono-posre.itp
sed -i "s/single-chain-10mono-posre/pp-single-chain-10mono-posre/" pp-single-chain-10mono.itp
We keep the unit cell dimensions unchanged for the nanofibril model ( 7 * 8 * 10.38 nm3) and fill the box with PP chains. To find the number of chains needed to fill the box with 0.85 g/cm3 density, we need to compute the empty space in the box.
We can calculate this information using the "gmx sasa" tool to compute the volume occupied by the cellulose nanofibril. First, we need to generate a temporary tpr file for input. The mdp parameter files (tarball file) for this tutorial are available as mdp-biocomp.tgz.
cp ~/Downloads/mdp-biocomp.tgz .
tar -xzf mdp-biocomp.tgz
gmx_mpi grompp -f em.mdp -p newtopol-mfc20.top -c mfc-20mono.gro -o sasa-mfc20
echo "Cellulose" | gmx_mpi sasa -f mfc-20mono.gro -s sasa-mfc20.tpr -o area-cellulose -tv volume-cellulose
Open the output file "volume-cellulose.xvg" using a text editor, and the last line of the file contains the volume information (see figure below).
The highlighted part in the last line marks the volume occupied by the cellulose nanofibril in the unit cell. We know the box volume is 581.28 nm3, so the volume available for the PP matrix will be 487.943 nm3. Now, we can compute the number of PP chains to be added to the box:
PP 10-monomer chain's mass = 422.796 amu or 702.0692136431 x 10⁻²⁴ g
Number of PP chains, X = (Density*Volume)/(Chain mass)
X = ~590 chains
To obtain the target density of 0.85 g/cm³, the composite box will require 590 PP chains of 10-monomer length.
To accommodate the addition of sufficient PP chains, we now generate a new box for the nanofibril system by extended the X and Y axes but keeping the Z-axis of the box the same as before to replicate an infinite (periodic) cellulose fibril and allow covalent links across the box boundary.
gmx_mpi editconf –f mfc-20mono.gro -o mfc-20mono-largebox.gro -c –box 18.0 18.0 10.38
Open the structure output "mfc-20mono-largebox.gro" and make sure that the fibril is centred and still has unchanged z-axis of the unit cell. Next, we insert the PP chains into the box:
gmx_mpi insert-molecules -f mfc-20mono-largebox.gro -ci pp-single-chain-10mono.gro -o biocomp-mfc20-pp10.gro -nmol 590 -scale 0.48 >& im-nmol590.log
#NOTE: The size of the new box has to be large enough to accommodate the required PP chains. In case the "insert-molecules" step does add sufficient molecules, increase the box size in X and Y dimensions and repeat the matrix insertion step until. Check the file im-nmol590.log to make sure enough chains are added for this model. If not, increase the box size and repeat the insert-molecules step.
Next, we create a new topology file using bash commands to include both cellulose and PP topologies.
sed "/Include water topology.*/i ; Include PolyPropylene chain topology\n#include \"pp-single-chain-10mono.itp\"\n " newtopol-mfc20.top >& newtopol.top
echo "Polypropylene_chain_X 590" >> newtopol.top
The structure model and topologies have been set up and we are ready to begin the MD simulations.
Energy Minimization
We use the same em.mdp parameter files as in the cellulose fibril tutorial for an optimization run. Now, run the same grompp and mdrun combination:
gmx_mpi grompp -f em.mdp -p newtopol.top -c biocomp-mfc20-pp10.gro -o em
gmx_mpi mdrun -deffnm em
The output structure from the minimization run will have some chains that are “broken” across the box. As before, reconstruct it using the trjconv tool:
echo "Cellulose System" | gmx trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro
Analyse the potential evolution of the system during the minimization run and evaluate if it's suitable to proceed to the equilibration step. (If not, repeat the minimization and try to check if there are any issues with the input structure).
Equilibration
By now it should be familiar that we begin with a short temperature equilibration with the V-rescale thermostat in the NVT ensemble.
gmx_mpi grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt
gmx_mpi mdrun -deffnm nvt -dlb yes
Upon successful completion, we again fix the “broken” output structure.
echo “Cellulose System” | gmx_mpi trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro
At this point, remember to check if the temperature is equilibrated to a reference value of 300 K. If not, then continue the NVT equilibration for another nanosecond.
Next, we implement the 2-step simulated annealing MD simulation (similar to the Polypropylene tutorial) to allow the box to compress and achieve a 0.85 g/cm3 density for the polypropylene matrix. The first run is a single annealing cycle over 500 ps (300 K → 450 K → 300 K). Before we begin the simulation, open the file npt-anneal-1.mdp and go through the thermostat, barostat and simulated annealing sections. You should notice that here we implement a semi-isotropic coupling for the barostat to allow scaling of box vectors only along X and Y axes.
gmx_mpi grompp -f npt-anneal-1.mdp -p newtopol.top -c nvt.gro -r nvt.gro -o npt-anneal-1
gmx_mpi mdrun -deffnm npt-anneal-1
(The "-dds" argument for mdrun allows the scaling of box vectors up to a ratio of 0.6)
Let's make the structure whole again using trjconv and visualize using PyMol.
echo “Cellulose System” | gmx_mpi trjconv -f npt-anneal-1.gro -s npt-anneal-1.tpr -pbc mol -center -o npt-anneal-1.gro
Open the structure file npt-anneal-1.gro and nvt.gro with PyMol. Notice the change in box vectors along the X and Y axes and the PP chains getting intertwined. Next, we perform some additional analysis to check the polymer density and cell size.
gmx_mpi energy -f npt-anneal-1.mdp -o npt-anneal-1.xvg
Choose the following options on the input prompt: Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, 0. Open the output xvg plot file and check the time evolution of each variable.
Next, start the second annealing MD run of 5ns time (npt-anneal-2.mdp) . During this run, the PP matrix will undergo 10 annealing cycles (300 K --> 600 K) with each repeat of 300 K → 600 K → 300 K occurring over 500 ps time. Here we again use separate temperature coupling for the cellulose fibril to keep it at the reference temperature of 300 K to avoid any structural deformation (see NVT section above). The mdp file parameters are similar to the previous step but now we add multiple heating/cooling cycles.
; SIMULATED ANNEALING
annealing = no single ; Type of annealing for each temperature group (no/single/periodic)
annealing-npoints = 2 21 ; Number of time points to use for specifying annealing in each group
annealing-time = 0 5000 0 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000 3250 3500 3750 4000 4250 4500 4750 5000
annealing-temp = 300 300 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300
Repeat the grompp and mdrun commands.
gmx_mpi grompp -f npt-anneal-2.mdp -p newtopol.top -r npt-anneal-1.gro -t npt-anneal-1.cpt -r npt-anneal-1.gro -o npt-anneal-2
gmx_mpi mdrun -deffnm npt-anneal-2 -dlb yes -dds 0.6
Again, we fix the broken output structure and evaluate the variables to check how the system fared against the annealing.
echo “Cellulose System” | gmx_mpi trjconv -f npt-anneal-2.gro -s npt-anneal-2.tpr -pbc mol -center -o npt-anneal-2.gro
gmx_mpi energy -f npt-anneal-2.edr -o npt-anneal-2.xvg
Choose the following options on the input prompt: Potential, Temperature, Volume, Density, 0. Open the output xvg plot file and check the time evolution of each variable.
The repeated annealing cycles will mould the PP chains to pack tightly around the fibril with steric interactions due to the hydrophilic-hydrophobic nature of the components. Here, it is important to remember that the densities stored in the 'npt-anneal-2.edr' file are for the whole composite system and not just the PP matrix. To measure the densities for PP and cellulose individually one can use the 'sasa' tool within gromacs. The left panel in the figure below shows the box sizes after each equilibration step (NVT: green, Anneal-1: yellow, Anneal-2: cyan) and the right panel shows the composite models before and after the annealing step. The structure are fitted to align the cellulose backbone atoms and depicts the tightly packed PP matrix in cyan.
We still need to verify that PP packing is stable and the box vectors will not deviate further under NPT conditions, so we perform another 10 ns equilibration run under NPT ensemble conditions temperature 300K and pressure 1.0 bar. This step also allows the PP matrix to equilibrate over a longer time at the reference temperature and stabilize the target densities. The parameter file for this run should already be in your current directory with the name 'npt-3.mdp'.
gmx_mpi grompp -f npt-3.mdp -p newtopol.top -c npt-anneal-2.gro -o npt-3
gmx_mpi mdrun -deffnm npt-3 -dlb yes
At the end of the run, use 'gmx energy' to analyze the box vectors to ensure stable volume and PP density. Also, remember to make the structure whole again using trjconv tool.