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.
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.
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.
