Preparation of initial configuration of polymers for MD simulation

It has been demonstrated that a single chain can be easily generated using the rotational isomeric state model. In a more complicated simulation with inclusion of the many-chain effect, assembly of many chains in a box and the corresponding equilibriation of the system can be a daunting task. But this can be easily achieved using free software, such as GROMACS and PACKMOL. In this tutorial, I will show the detailed procedures for preparing a system of polyethylene chains with 20 carbon atoms. First and foremost, we need to prepare a protein database bank file (pdb) of a single chain. Such a pdb file can be downloaded in the following link. (initial.pdb) Such conformation was generated by the rotational isomeric state model, which accounts for the interdependence of the trans and gauche states.

Secondly, we then use PACKMOL software to pack 40 chains in a box with dimension $3.0\times 3.0\times 3.0$ nm$^3$.

# # PE system # # All the atoms from diferent molecules will be separated at least 2.0 # Angstroms tolerance 2.0 # The file type of input and output files is PDB filetype pdb # The name of the output file output initial-40.pdb structure initial.pdb number 40 inside box 0. 0. 0. 30. 30. 30. radius 2 end structure

Run the following command in the terminal:

packmol < pack_pe.inp

After that, we have to properly convert the initial-40.pdb file to a pdb file that can be understood by GROMACS, which can be achieved using the following Python 3.5 program (pack2pdbpe). A pdb file, named initial2.pdb, will be generated by that program.

Thirdly, we need to let GROMACS know the atom type and the definitions of PEM, PEE as well as PEN in the pdb file. Therefore, we need to add in the following to the aminoacids.rtp file.

[ PEE ] [ atoms ] C1 opls_966 0.000 1 [ bonds ] -C1 C1 C1 +C1 [ PEN ] [ atoms ] C1 opls_967 0.000 1 [ bonds ] -C1 C1 [ PEM ] [ atoms ] C1 opls_967 0.000 1 [ bonds ] C1 +C1

Fourthly, the following command can be run:

gmx pdb2gmx -f initial2.pdb

This will generate several more files, such as conf.gro, topol.top, etc., which are necessary in the energy minimization and the actual molecular dynamics simulation. Now, we can perform energy minimization with the file em.mdp. But first, let us generate a box for the system.

gmx editconf -f conf.gro -box 3.3 -o conf2.gro

Note that we choose a box with a dimension slightly larger than that used in the PACKMOL. This is because PACKMOL software does not include periodic boundary conditions. With these, we can generate a .tpr file for the energy minimization.

gmx grompp -f em.mdp -p topol.top -c conf2.gro -o em.tpr

The em.mdp file is as follows:

title = fws cpp = /usr/bin/cpp define = -DFLEX_SPC constraints = none integrator = steep dt = 0.1 ; ps ! nsteps = 5000 nstlist = 10 nstenergy = 500 ns_type = grid r_list = 0.7 ; ; Energy minimizing stuff ; emtol = 100.0 emstep = 0.01

Now, run the energy minimization with the following command.

gmx mdrun -v -deffnm em

This is then followed by the compression of the system in a NPT simulation with the pressure coupling using Berendsen barostat. The md.mdp file can be seen here.

; RUN CONTROL PARAMETERS integrator = md tinit = 0 dt = 0.001 ; 1 femtosecond time step for integration nsteps = 10000000 ; Make it 10ns to produce the traj. ; OUTPUT CONTROL OPTIONS nstxout = 500000 nstvout = 500000 ;nstfout = 100 ;nstlog = 100 ;nstenergy = 100 ;nstxtcout = 100 xtc_precision = 100 xtc-grps = System energygrps = System ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns-type = grid pbc = xyz rlist = 0.9 ; OPTIONS FOR ELECTROSTATICS AND VDW ;coulombtype = PME rcoulomb = 1.0 fourierspacing = 0.12 vdw-type = Cut-off rvdw = 1.0 ; Temperature coupling tcoupl = nose-hoover tc-grps = system tau_t = 0.2 ref_t = 450 ; Pressure coupling Pcoupl = Berendsen Pcoupltype = Isotropic tau_p = 0.5 compressibility = 4.5e-5 4.5e-5 4.5e-5 ref_p = 1.01325 ;bar ; OPTIONS FOR BONDS constraints = none

Similar to what has been done for the energy minimization, we run the follow command to generate a .tpr file, and then run the simulation.

gmx grompp -f md.mdp -p topol.top -c em.gro -o md.tpr -maxwarn 3 gmx mdrun -v -deffnm md

After such equilibration, we can obtain a reasonable initial configuration for simulation in NVT ensemble.

Figure 1: Initial configuration of many-chain system for NVT simulation.

Go back; Go back to Main Page