Polypropylene matrix MD tutorial

(This tutorial is related to work published at https://doi.org/10.3390/nano12193379).

In this section, we will generate a polypropylene (PP) matrix and perform preliminary calculations to analyze the characteristics of the box. We start by building the structure of a 10-monomer propylene chain using the CHARMM-GUI interface. The procedure to build a polymer model are described in a video demo. Here are the input choices to be selected to build a 10-monomer propylene chain using the Polymer Builder tool: 

Polymer Builder --> (System Type) Single Polymer Chain --> (Building block) Olifens --> Polypropylene --> isotactic [10] 

After selecting the above options, click on the "Build Polymer Chains" button. Once the job is completed in few minutes, download the PDB structure file named “p1_raw.pdb” by clicking on the filename. To obtain a complete overview of the methodology for the polymer structure generation in CHARMM-GUI, read a recent article by Choi YK et.al. 

Alternatively, you can also the download the structure model for a single chain of PP with 10-monomers.

Next, we use a set of bash commands to modify the residue and atom names that are compatible with the RTP description of PP in Gromacs. It is necessary to note that in this PP chain, the terminal residues would have a different number of atoms in comparison to the internal residues. To account for this in the molecule description in Gromacs, we use a three-residue approach implemented for polymer MD simulations by Wildman, J et.al. 

For the RTP file description in force fields, we use three residues to emulate a starting residue (SPP), 8x internal residues (IPP), and an ending residue (EPP) of the chain. At this point, you can open the structure file using PyMol to internal and terminal residues in the chain. Next, we edit the atom and residue names in the newly built PP structure file to match the description in the RTP file of force fields.

mkdir PP-10mono 

cp ~/Downloads/p1_raw.pdb PP-10mono/pp-10mono_raw.pdb 

cd PP-10mono 

sed 's/PROR /IPP X/; s/IPP X   1/SPP X   1/; s/IPP X  10/EPP X  10/; /IPP/s/H13 /H12 /; /EPP/s/H13 /H12 / ; s/TER.*/TER/' pp-10mono_raw.pdb >& pp-10mono.pdb 

The propylene molecule is not a standard amino acid, the force field parameters for Amber14sb were derived for the polymer chain. The single-chain structure was optimized using the DFT level of theory (HF/6-31G*). The procedure has been adapted based on developer comments in the AMBER forum for parameterizing polymer chain residues when using the RESP method (https://doi.org/10.1002/jcc.540150705).

 Next, we start with generating the topology and Gromacs structure file for the PP chain using the Amber14 force fields. 

echo " 1" | gmx_mpi pdb2gmx -water tip3p -f pp-10mono.pdb -o single-chain-10mono.gro -p topol-pp.top -i single-chain-10mono-posre.itp 

If you are curious about the force field parameter files for the PP chain, go to the amber14sb_OL15.ff subdirectory under the force-fields directory and check the files with the “pp” prefix. 

In the next step, we want to create a box and fill it with PP chains to form the matrix. Instead of generating a topology file for the chain individually, we can use the single-chain description as a repeating unit like adding water to the cellulose fibril box in the previous section. For this we first convert the “top” file into an “itp” format topology file to use as a base for describing the multiple chains in the box and create a new topology file “newtopol.top”. 

sed -n '/moleculetype/,/endif/p' topol-pp.top | sed 's/Composite_chain_X/Composite        /' >& single-chain-10mono.itp 

sed '/moleculetype/,/endif/d' topol-pp.top | sed "/Include water topology.*/i ; Include PolyPropylene chain topology\n#include \"single-chain-10mono.itp\"\n " >& newtopol.top 

Now we can use the Gromacs tool insert-molecules tool to add the chains randomly into a box of fixed size. However, the random or unordered insertion of the wooden log like chains into a box means we cannot at once obtain the target density without avoiding steric clashes. The solution to this problem is to first compute the number of chains we would need in an 8*8*8 nm3 box. Next, we add the computed chain count to a much larger box (12*12*12 nm3) and later use a simulated annealing approach to allow the box to compress to target volume/density. Below is the formulation on how to compute the number of 10-monomer polypropylene chains we would need to achieve ~0.85 g/cm3 in an 8*8*8 nm3 box. 

Atomic mass of 10-monomer chain: 422.796 amu or 702.0692136431 x 10⁻²⁴ g 
Unit cell volume: 512 nm³ or 5.12 * 10-19 cm³
Target chain count X (10-monomer chain): (Density*Volume)/(Chain mass) 
For our system, X: (0.85 g/cm³ * 5.12* 10-19 (cm³))/ (702.0692136431 * 10⁻²⁴ g) = ~620 chains 

Use the following Gromacs command to create the large box with ~620 polymer chains: 

gmx_mpi insert-molecules -ci single-chain-10mono.gro -box 12 12 12 -nmol 620 -scale 0.48 -o box-large.gro >& inmol-620chains.log  

Note that the default vdw radii value of “0.57” is meant for a density of 1 g/cm3 for water, but here we will use a “0.48” radii value to obtain ~0.85 g/cm3 density. 

Important: In practice, the size of the large box is decided by a trial-and-error approach until we find the box size that fits the target polymer chain count. This approach is based on the study by Lemkul. JA et.al.  

Now we edit the [ molecules ] section of the topology file “newtopol.top” to add that we have 620 PP chains in the system. 

[ molecules ] 
; Compound        #mols 
Polypropylene_chain_X  620 

Next, run a short energy minimization with 1000 kJ/mol threshold. We can use the same em.mdp file but need to comment out the line for periodic_molecules since have finite chain length here. 

; periodic_molecules = yes

Repeat the same grompp and mdrun combination now: 

gmx_mpi grompp -f em.mdp -p newtopol.top -c box-large.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 “Polypropylene System” | gmx_mpi trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro 

Now analyze the potential energy curve of the minimization run using the same “gmx energy” tool in Gromacs. 

Now start the 1 ns temperature equilibration run in the same manner as before using the parameter file nvt.mdp but with few modifications in the parameters: 

;periodic_molecues = yes 
; Temperature coupling 
tcoupl                   = V-rescale    ; V-rescale Thermostat 
tc-grps                  = System 
tau-t                    = 0.1          ; coupling time constant [ps] 
ref-t                    = 300          ; reference temperature [K] 

Run the following Gromacs commands to start the temperature equilibration run.

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 “Polypropylene System” | gmx_mpi trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro 

Additionally, use the energy tool to check the temperature evolution and visualize the output structure in the file nvt.gro. You will now see that PP chains start to bend and form intermolecular interactions. 

Now that the temperature is up to reference value, we equilibrate the pressure and include annealing parameters to take the temperature of the system from 300 K to 450 K and back to 300 K. We use an isotropic coupling here to allow uniform scaling of the three-box vectors, otherwise, we may end up with a non-cubic (cuboid) box. The parameters modified are shown below: 

; Pressure coupling 
pcoupl                  = Berendsen     ; Pressure coupling on in NPT 
pcoupltype              = isotropic     ; uniform scaling of box vectors 
tau_p                   = 1.0           ; time constant, in ps 
ref_p                   = 1.0           ; reference pressure, in bar
compressibility         = 4.5e-5        ; isothermal compressibility of water, bar^-1

; SIMULATED ANNEALING 
annealing                = single                   ; Type of annealing for each temperature group (no/single/periodic) 
annealing-npoints        = 6                        ; Number of time points to use for specifying annealing in each group 
annealing-time           = 0 100 200 300 400 500    ; time in ps for annealing points for each group  
annealing-temp           = 300 300 375 450 375 300  ; temperature at the annealing point. 

You can download the parameter file npt-anneal-1.mdp. 

gmx_mpi grompp -f npt-anneal-1.mdp -p newtopol.top -c nvt.gro -r nvt.gro -t nvt.cpt -o npt-anneal-1 

gmx_mpi mdrun -deffnm npt-anneal-1 -dlb yes -dds 0.6 

Notice the additional specification of "dds" value for the mdrun run input. The value dictates the reciprocal of a fraction by which the cell size can increase/decrease during the run while preserving the domain decomposition of the system cell. Without this parameter, the annealing run is very likely to crash due to the compression of the box size. 

Before looking at the output structure, make sure to repair the “broken” structure using trjconv 

echo “Polypropylene 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. The matrix cell shrunk significantly after the annealing NPT run and the chains are fully intertwined as expected from an actual PP matrix. Next, we perform some additional analysis to check the polymer density and cell size. 

 gmx_mpi energy -f npt-anneal-1.edr -o npt-anneal-1.xvg 

Choose the following options on the input prompt: Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, Density, 0. Open the output xvg plot file and check the time evolution of each variable. 

Next, start another annealing MD of 5 ns time with the new box. Now we improve the polymer equilibration by heating the system from 300K --> 600 K about 5-10 times with each repeat of 300K-->600K-->300K occurring over 500ps time. In some cases, that system cell may not compress or might still be far from the target density, hence the second annealing cycle. The mdp file parameters are similar to the previous step but now we add multiple heating/cooling cycles. 

; SIMULATED ANNEALING 
annealing                = single       ; Type of annealing for each temperature group (no/single/periodic) 
annealing-npoints        = 21           ; Number of time points to use for specifying annealing in each group 
annealing-time           = 0 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000 3250 3500 3750 4000 4250 4500 4750 5000 
annealing-temp           = 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 600 300 

Download the parameter file npt-anneal-2.mdp and repeat the grompp and mdrun commands.

gmx_mpi grompp -f npt-anneal-2.mdp -p newtopol.top -c 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 “Polypropylene 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, Pressure, Box-X, Box-Y, Box-Z, Volume, Density, 0. Open the output xvg plot file and check the time evolution of each variable. Below is an example plot from a similar run for a 7*7*7 nm3 box of PP. 

 Next, we run an additional 10 ns NPT equilibration (without annealing and position restraints) to stabilize the new cell. With this final equilibration run, the aim is to stabilize the cell size and atomic interactions. The parameter file for this step is npt-3.mdp. 

gmx_mpi grompp -f npt-3.mdp -p newtopol.top -c npt-anneal-2.gro -t npt-anneal-2.cpt -r npt-anneal-2.gro -o npt-3 

gmx_mpi mdrun -deffnm npt-3 -dlb yes 

In case the simulation crashes then restart the run with a smaller time-step of 1-fs to avoid any system blow-up issues. When the run ends successfully, run the trjconv and energy commands with appropriate arguments to make the output structure (npt-3.gro) whole again and to analyze the stability of box vectors, volume and PP density.