I’ve been trying to perform steered molecular dynamics simulations and have a couple of questions.
I do not have and have not had good support from my supervisor, so classify me as having very basic knowledge. I would really appreciate some help or a pointing towards a good page/tutorial/video or something else.
I’m using pulling groups (as a teacher called it) and umbrella sampling to pull an ion from a pair of amino acid side chains.
My issue has been that the ion in the frame used for the SMD is between the side chains. But in the simulation it stays glued to the oxygens in the COO- end of one of the side chains.
I heard from a person who really didn’t have time to help me (understandably) very short to look into PLUMED. Does someone know how PLUMED might be able to help?
The other thing is that I don’t know how the force-time or displacement-time plots should look. I get the below but I don’t know either how to interpret them or if they are correct looking.
I’m not sure if I’m giving you what you were after, so just ask again if you need something more/else.
The mdrun md.mdp-file content is below.
Okay so it is not like PLUMED can specify pulling groups in another way, or does the equilibrium, nvt and npt runs which I do before the mdrun, in a different way where I end up with the pulling groups being like in my starting frame - i.e. a pair of amino acid residue with the ion between the COO- ends of their side chains?
What is actually the difference to PLUMED from using pulling groups the way I am?
md.mdp:
title = beta-galactosidase amber99sb-ildn tip3p
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 1000000
dt = 0.002 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; 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 is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; 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 is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in 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 ; Velocity generation is off
; COM pulling
pull = yes
pull-print-ref-value = yes
pull-ncoords = 2
pull-ngroups = 3
pull-group1-name = a_3999-4000
pull-group2-name = a_18121-18122
pull-group3-name = a_32615
pull-coord1-groups = 1 3
pull-coord2-groups = 2 3
pull-coord1-type = umbrella
pull-coord2-type = umbrella
pull-coord1-geometry = distance
pull-coord2-geometry = distance
pull-coord1-init = 0.0
pull-coord2-init = 0.0
pull-coord1-start = yes
pull-coord2-start = yes
pull-coord1-rate = 0.00075 ; nm/ps
pull-coord2-rate = 0.00075 ; nm/ps
pull-coord1-k = 3000 ; kj/mol/nm^2
pull-coord2-k = 3000 ; kj/mol/nm^2
Additional information:
When I view the trajectory in VMD it does look like the distance increases to the magnesium ion due to the side chains of the amino acids bending. It looks like the magnesium is attached to one or two of the oxygens (depending on if I use 4 groups for the oxygens or 2, and depending on the force constant I use) and that either this “combo” is being pulled from the other oxygen pair, or the opposite.
What I can’t understand is how the force and distance isn’t proportional to each other. Shouldn’t they follow the same pattern (one goes up, so does the other) in umbrella sampling?
For the bending of aminoacids: If I instead had chosen an atom such as the first carbon of the side chain, would the magnesium still be feeling the charge of the COO- group at the top of the side chain? And would the magnesium simply be somewhere around one of the first carbons, instead of at the COO- group’s charge?
The force goes up because the actual distance does not follow the reference distance. Something is hindering the distance from increasing. I can’t say what that is.
Note that the distance you show is very large: you are 1 nm away from the initial distance. I would expect the ion to be in solution and not be hindered at all. Or is it stuck to some other part of the molecule?
