[gmx-users] Gromacs - Dimer Simulation problems

[Michael Carter]

Hi,

I have been running a short MD simulation (100ns) on a dimer protein. The two monomers are identical, after looking into your suggestions on http://www.gromacs.org/Documentation/How-tos/Multiple_Chains I  have set up the system with waters and ions. The topol.top to look like this:

;
;       File 'topol.top' was generated
;       By user: michaelc (1128283773)
;       On host: 104844CTHDT.local
;       At date: Tue Jul  1 12:30:36 2014
;
;       This is a standalone topology file
;
;       It was generated using program:
;       pdb2gmx - VERSION 4.6.5
;
;       Command line was:
;       /usr/local/gromacs/bin/pdb2gmx -f BCL6_584_pseudo_apo.pdb -o BCL6_584_pseudo_apo.gro -water tip3p -ignh
;
;       Force field was read from the standard Gromacs share directory.
;

; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"

; Include chain topologies
#include "topol_Protein_chain_A.itp"
#include "topol_Protein_chain_B.itp"

; Include water topology
#include "amber99sb.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include topology for ions
#include "amber99sb.ff/ions.itp"

[ system ]
; Name
BCL6_584_pseudo_apo in water

[ molecules ]
; Compound        #mols
Protein_chain_A     1
Protein_chain_B     1
SOL         40832
NA               76
CL               79

I have then run a short minimisation and equilibration which ran successfully. I then ran a production MD run of 100ns. The md.mdp looks like this:

title       = Protein-ligand complex MD simulation
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 50000000  ; 2 * 50000000 = 100000 ps (100 ns)
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 0         ; suppress .trr output
nstvout     = 0         ; suppress .trr output
nstenergy   = 1000      ; save energies every 2 ps
nstlog      = 1000      ; update log file every 2 ps
nstxtcout   = 1000      ; write .xtc trajectory every 2 ps
energygrps  = System
; Bond parameters
continuation    = yes           ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
ns_type     = grid      ; search neighboring grid cells
nstlist     = 5         ; 10 fs
rlist       = 0.9       ; short-range neighborlist cutoff (in nm)
rcoulomb    = 0.9       ; short-range electrostatic cutoff (in nm)
rvdw        = 1.4       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = Protein Water_and_ions    ; two coupling groups - more accurate
tau_t       = 0.1   0.1                     ; time constant, in ps
ref_t       = 300   300                     ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl      = Parrinello-Rahman             ; pressure coupling is on for NPT
pcoupltype  = isotropic                     ; uniform scaling of box vectors
tau_p       = 2.0                           ; time constant, in ps
ref_p       = 1.0                           ; reference pressure, in bar
compressibility = 4.5e-5                    ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; assign velocities from Maxwell distribution

After this simulation finished I obtained the resultant .xtc and .tpr files. The system was then corrected to remove any problems of drifting cause by PBC and the system was centered on the protein.  This was done using the following commands:


trjconv -s md_0_100.tpr -f md_0_100.xtc -o protein.xtc
tpbconv -s md0_100.tpr -o protein.tpr

Choosing the group 1 (Protein).

trjconv -s protein.tpr -f protein.xtc -o protein.gro -dump 0

Choosing the group 1 (Protein).

trjconv -s protein.tpr -f protein.xtc -o protein_fit.xtc -fit rot+trans

Chossing the group 4 (Backbone) and the system 1 (Protein) to write to disk.

After analysing the output it appears that after around 50ns something very strange happens and my simulation see a large jump in RMSD to over 5A. This continues to happen for the last 40ns of the simulation. Upon visualisation it would appear that the dimer dissociates into two monomers. There is no real reason for this and it appears to happen randomly for the last 50 ns of the simulation.

Have you come across this problem before? Do you think that I should apply some constraints to stop this from happening? Or could it be that I need to adjust my system setup to remove this issue?

Any help on this matter would be greatly appreciated.

All the best,
Michael Carter

The Institute of Cancer Research: Royal Cancer Hospital, a charitable Company Limited by Guarantee, Registered in England under Company No. 534147 with its Registered Office at 123 Old Brompton Road, London SW7 3RP.

This e-mail message is confidential and for use by the addressee only.  If the message is received by anyone other than the addressee, please return the message to the sender by replying to it and then delete the message from your computer and network.