[gmx-users] Effect of refcoord_scaling

[Matthias Ernst]

Hi Chris,

thank you for your answer.
Let me comment on some of your hints.
> refcoord_scaling is only required when you are also using positions restraints.
> Therefore we need to know what exactly you are doing with position restraints
> in order to provide the most useful advice.
Yup, that's right, I forgot to mention this (sorry). I am running 
simulations of a complex of a protein with a DNA double helix. In test 
simulations, I found that the DNA distorts almost immediately if it's 
free (it looks like the double helix will un-wind and bend, resulting in 
a lot of mdrun warnings and finally abort because of large movement, 
even at 0.5fs timestep) and this movement affects the DNA-protein 
interactions.
To avoid the distortion, I thought I can apply position restraints on 
the DNA to keep it in place (more or less) to get a better picture of 
the interaction. Is there another way to do this?

And what came to my mind when I considered this, if I apply position 
restraints to a molecule to kind of "fix it" in a NpT simulation, should 
I include or exclude it in the tc_grps (or maybe include at T=0)?
> Nevertheless, you ran 2 simulations and got different results. It is not prudent to
> assign the difference to refcoord_scaling at this point. To test this yourself, please
> repeat each simulation (ideally at least 3 simulations for each case with and
> without refcoord_scaling).
That sounds reasonable and I think I ought to do this. Unfortunatly the 
simulations are quite expensive and I have a (quite hard) deadline next 
week when I have to submit my diploma thesis. So, I think I will not be 
able to repeat these simulations before. Even if I start them now, they 
will not be finished (even if they start soon which, too, is unlikely). 
Still, I can do this afterwards.
> I am not sure what happens with pressure coupling but using refcoord_scaling=no
> (the default). The manual says:
> "Note that with this option the virial and pressure will depend on the absolute
> positions of the reference coordinates."
> I interpret this to mean that you will get the wrong pressure, and my hunch is that
> this would not significantly affect the stability of a DNA-protein complex, but you'll need
> to test that out yourself.
This is exactly the point why I wanted to ask if somebody has 
experiences with this issue and can tell us what this combination of 
parameters can cause in a simulation.
> A final note is that you should be sure to use the exact same conformation to start your
> runs both with and without refcoord_scaling=com. Either start with this conformation and
> redo the minimization, solvation, etc for each replica or pick one of your minimized initial
> conformations to start all of your production runs. This is important so that you avoid
> the situation in which some stochastic event in your system setup (pre production runs)
> actually lead to the difference.
Okay. What I did was to start with the same structure and then apply a 
several-step protocol similar to the tutorials to it: EM in vacuo, add 
solvent and EM, add ions (only to neutralize it, no excess salt 
concentration) and EM, MD with entire system position-restrained, MD in 
NVT-ensemble for equilibration of temperature, MD in NpT-ensemble for 
equilibration of pressure and then, finally, production MD. This whole 
protocol was  carried out for both of my simulation, so the initial 
positions of the ions are quite different and maybe this plays a role, 
too. Apart from this (and the velocities of the particles), the setup is 
identical.
The proper step to "jump in" when repeating the simulations seems to be 
before NpT-equilibration.

Please, if you see some obvious (or not so obvious) mistake in what I'm 
doing, tell me. One question that also could not be resolved after 
reading the tutorials was if it is good/necessary or rather a bad idea 
to continue with the old velocities in the equilibration steps. Some 
tutorials do, others don't...

Sorry if there are some rather simple questions, but unfortunately I 
don't have a supervisor who knows GROMACS and who could tell me what to 
do or answer my questions. In addition, I did not have much time to get 
used to all this as in the beginning, the project was meant to use MC 
simulations instead of MD what took me a rather long time to implement 
and to find out that this does not work well.

But still, I have some results. I want to understand them well and make 
sure next time, I will do less mistakes...

Thanks for your help,
Matthias

> Chris.
>
> -- original message --
>
> I am currently working on Protein-DNA-complexes. They should be
> simulated in NPT-ensemble.
> I did the same simulation including previous minimization steps (in
> vacuo, with solvent, with solvent and ions) and equilibration (system
> position restrained, with theromstate, with barostate) twice with one
> slight difference: in the second case, there was a GROMPP warning that
> NPT (Berendsen-barostate) needs refcoord_scaling to avoid artifacts,
> therefore I added "refcoord_scaling=com" to my mdp file.
> The systems showed significantly different behaviour. In the first run
> (without refcoord_scaling), the protein-DNA-complex was unstable and
> some of the contacts between them broke. In the second run, the complex
> remained stable.
> As I do not have much experience with explicit solvent and ions MD
> simulations, wondered if this difference can be caused by the lack of
> reffcoord_scaling command.
> The other guess would be that this comes due to an ion that drifts in
> between the DNA and and the protein and therefore causes the distortion
> of the protein.
>
> Which do you think would be more likely? And which types of artifacts
> can be caused by lack of refcoord_scaling and can they be seen or
> detected easily?
>
> Thank you very much for your help,
> Matthias Ernst