GROMACS version: 2020
GROMACS modification: No
I’m conducting a coarse-grained molecular dynamics simulation where I’m pulling my peptide through to the centre of mass of my bilayer system. The distance between the protein’s COM and membrane is ~ 3nm, the time is 50ns, so my pulling rate is about 0.0007 nm/ps with a force constant of 1,000. When I plot my force evolving with time it begins from 0 then tends to negative kJ/mol/nm, could anyone help with why this is? I feel like my system is unrealistic but my pulling force feels pretty reasonable.
I should add, when the peptide first reaches the membrane I get my maximum on the PMF free energy which I assume is my transition state, I’m guessing my protein has just gone further than the membrane’s COM so the force is now acting in the opposite direction to restore it towards the COM but not sure how to avoid this situation?
It is very hard to conclude anything from a non-equilibrium simulation. You should umbrella potentials at fixed locations, or AWH which is more efficient and easier to set up.
I’m using distance between the COM of the protein and the membrane. Is this the starting step for the umbrella potentials at fixed locations or is it something I would put in the input .mdp file?
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = Protein
pull_group2_name = Membrane
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N N Y
pull_coord1_groups = 1 2
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = -0.00007 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
