[gmx-users] Effect of refcoord_scaling

[Matthias Ernst]

Hi Erik,

I guess you want to see the production mdp file. Here it is:


; VARIOUS PREPROCESSING OPTIONS
title                    = Production Simulation
cpp                      = /lib/cpp
define                   = -DPOSRES_DNA

; RUN CONTROL PARAMETERS
integrator               = md
tinit                    = 0       ; Starting time
dt                       = 0.0005 ; 2 femtosecond time step for integration
nsteps                   = 100000000

; OUTPUT CONTROL OPTIONS
nstxout                  = 50000 ; Writing full precision coordinates 
every nanosecond
nstvout                  = 50000 ; Writing velocities every nanosecond
nstfout                  = 0     ; Not writing forces
nstlog                   = 5000  ; Writing to the log file every step
nstenergy                = 5000  ; Writing out energy information every step
nstxtcout                = 5000  ; Writing coordinates every 5 ps
energygrps               = Protein DNA water_and_ions

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 5
ns-type                  = Grid
pbc                      = xyz
rlist                    = 0.9

; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype              = Reaction-Field
rcoulomb                 = 1.4
epsilon_rf               = 78
epsilon_r                = 1
vdw-type                 = Cut-off
rvdw                     = 1.4

; Temperature coupling
Tcoupl                   = Berendsen
tc-grps                  = Protein  DNA   water_and_ions
tau_t                    = 0.1      0.1      0.1
ref_t                    = 200      200   200
; Pressure coupling
Pcoupl                   = Berendsen
Pcoupltype               = Isotropic
tau_p                    = 1.0
compressibility          = 4.5e-5
ref_p                    = 1.0

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel                  = no

; OPTIONS FOR BONDS
constraints              = all-bonds
constraint-algorithm     = Lincs
unconstrained-start      = yes
lincs-order              = 4
lincs-iter               = 1
lincs-warnangle          = 30


Thanks for your help,
Matthias


Am 06.09.2012 09:47, schrieb Erik Marklund:
> Hi,
>
> Although your preparations seem sensible and a good starting point for stable simulations, I suspect that the source of the error is in your mdp-file or in how you constructed your simulation box. I've simulated complexes like yours for hundreds of nanoseconds without problems. Could you provide your mdp-file for us to look at?
>
> Best,
>
> Erik
>
> 6 sep 2012 kl. 00.31 skrev 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
>> -- 
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-request at gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> -----------------------------------------------
> Erik Marklund, PhD
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,    75124 Uppsala, Sweden
> phone:    +46 18 471 6688        fax: +46 18 511 755
> erikm at xray.bmc.uu.se
> http://www2.icm.uu.se/molbio/elflab/index.html
>