MD simulation of system of benzene rings

Alternatively, one may not even need to use PACKMOL for initial configuration preparation, but only by using the insert-molecules commanded provided by GROMACS-5.1.4. The procedures are extremely straightforward and easy. Insert a benzene fragment in Avogadro (Click Build>Insert>Fragment>aromatics>benzene.cml). Then, click Build>Remove Hydrogens>Add Hydrogens. With this step, the coordinates of the carbon and hydrogen atoms in the PDB file will be more easily managed as one will see in the forthcoming discussion.

The PDB file is explicitly shown as follows:

COMPND UNNAMED AUTHOR GENERATED BY OPEN BABEL 2.3.2 ATOM 1 C LIG 1 0.682 -0.092 1.209 1.00 0.00 C ATOM 2 C LIG 1 -0.707 -0.035 1.197 1.00 0.00 C ATOM 3 C LIG 1 -1.390 0.057 -0.011 1.00 0.00 C ATOM 4 C LIG 1 -0.682 0.092 -1.209 1.00 0.00 C ATOM 5 C LIG 1 0.707 0.035 -1.197 1.00 0.00 C ATOM 6 C LIG 1 1.390 -0.057 0.011 1.00 0.00 C HETATM 7 H LIG 1 1.188 -0.160 2.106 1.00 0.00 H HETATM 8 H LIG 1 -1.232 -0.061 2.085 1.00 0.00 H HETATM 9 H LIG 1 -2.421 0.099 -0.019 1.00 0.00 H HETATM 10 H LIG 1 -1.188 0.160 -2.106 1.00 0.00 H HETATM 11 H LIG 1 1.232 0.061 -2.085 1.00 0.00 H HETATM 12 H LIG 1 2.421 -0.099 0.019 1.00 0.00 H CONECT 1 2 6 7 CONECT 2 1 3 8 CONECT 3 2 4 9 CONECT 4 3 5 10 CONECT 5 4 6 11 CONECT 6 1 5 12 CONECT 7 1 CONECT 8 2 CONECT 9 3 CONECT 10 4 CONECT 11 5 CONECT 12 6 MASTER 0 0 0 0 0 0 0 0 12 0 12 0 END

Note that the coordinates of the carbon atoms and hydrogen atoms are organized in a way that the first 6 rows include the coordinates of the carbon atoms, followed by that of the hydrogen atoms. Nonetheless, a PDB file in such format cannot be read by GROMACS. Therefore, a Python program was written to regenerate a PDB file that is readable by GROMACS. (the program can be downloaded here: pack2pdb) The final PDB file for simulation is shown as follows:

ATOM 1 C1 BEC 1 18.04 16.5 10.4 1.00 0.00 ATOM 2 C2 BEC 1 19.01 17.41 9.99 1.00 0.00 ATOM 3 C3 BEC 1 18.64 18.68 9.59 1.00 0.00 ATOM 4 C4 BEC 1 17.3 19.06 9.6 1.00 0.00 ATOM 5 C5 BEC 1 16.33 18.16 10.02 1.00 0.00 ATOM 6 C6 BEC 1 16.7 16.88 10.42 1.00 0.00 ATOM 7 H1 BEC 1 18.32 15.55 10.7 1.00 0.00 ATOM 8 H2 BEC 1 20.0 17.13 9.97 1.00 0.00 ATOM 9 H3 BEC 1 19.36 19.36 9.28 1.00 0.00 ATOM 10 H4 BEC 1 17.03 20.01 9.31 1.00 0.00 ATOM 11 H5 BEC 1 15.34 18.43 10.03 1.00 0.00 ATOM 12 H6 BEC 1 15.98 16.2 10.73 1.00 0.00

Before importing the PDB file into GROMACS, we edit the forcefield file, which is named aminoacids.rtp in the folder oplsaa.ff:

[ BEC ] [ atoms ] C1 opls_145 0.000 1 C2 opls_145 0.000 1 C3 opls_145 0.000 1 C4 opls_145 0.000 1 C5 opls_145 0.000 1 C6 opls_145 0.000 1 H1 opls_146 0.000 1 H2 opls_146 0.000 1 H3 opls_146 0.000 1 H4 opls_146 0.000 1 H5 opls_146 0.000 1 H6 opls_146 0.000 1 [ bonds ] C1 C2 C2 C3 C3 C4 C4 C5 C5 C6 C6 C1 C1 H1 C2 H2 C3 H3 C4 H4 C5 H5 C6 H6 [ impropers ] C1 C2 C3 C4 improper_Z_CA_X_Y C6 C1 C2 C3 improper_Z_CA_X_Y C5 C6 C1 C2 improper_Z_CA_X_Y C4 C5 C6 C1 improper_Z_CA_X_Y C3 C4 C5 C6 improper_Z_CA_X_Y C2 C3 C4 C5 improper_Z_CA_X_Y

Now, run the command:

gmx pdb2gmx -f initial2.pdb

This will generate several files, that are necessary for MD simulation. Let us generate a box with 40 benzene rings:

gmx insert-molecules -ci conf.gro -box 2.2 2.2 2.2 -nmol 40 -o conf2.gro

The energy of the system is then minimized (More information of the energy minimization can be found on Preparation of initial configuration of polymers for MD simulation). However, a system generated in such a way is not as perfect as that in PACKMOL. Therefore, we firstly need to run a NVT simulation before compressing the system to the correct density in the NPT simulation. The md.mdp file is as follows:

; RUN CONTROL PARAMETERS integrator = md tinit = 0 dt = 0.001 ; 1 femtosecond time step for integration nsteps = 500000 ; Make it 10ns to produce the traj. ; OUTPUT CONTROL OPTIONS nstxout = 1000 nstvout = 1000 ;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 = 273 ; 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 ; GENERATE VELOCITIES FOR STARTUP RUN ;gen_vel = yes ; ;gen_temp = 450 ; ;gen_seed = 173529 ; ; OPTIONS FOR BONDS constraints = none

We can then run the following command:

gmx grompp -f md.mdp -p topol.top -c em.gro -o md.tpr -maxwarn 1

Followed by:

gmx mdrun -v -deffnm md

As shown in this md.mdp file, the pressure coupling section is commented with the sign ; such that pressure coupling is disabled. But when we perform NPT simulation in the future, this section has to be uncommented. Surprisingly, the system can be equilibriated within a short period time (i.e. 500 ps of NVT simulation and 500 ps of NPT simulation).

Figure 1: Final configuration of 40 benzene rings after 500 ps of NVT and 500 ps of NPT simulations.

If you think this demonstration is unclear, please find the following video.

Go back; Go back to Main Page