[gmx-users] Effect of refcoord_scaling

Christopher Neale chris.neale at mail.utoronto.ca
Thu Sep 6 18:15:17 CEST 2012

What you should do depends on your question. If you question is represented by your original post and by your subject (what is the effect of refcoord_scaling?) then you need to repeat both conditions a few times. If you just want to do it correctly, then use refcoord_scaling=com (note that this is not a guarantee of success if you have setup your system incorrectly or if your parameters lead to an inherently unstable association).

If your system is unstable at 300 K, then I do not agree that this makes it acceptable to run at 200 K. This may just be hiding some underlying problem with your setup. For example, if you simulate at 200 K for a while and then move to 300 K, does that lead to a stable system? If so, then the problem was with your initial conformation. What has been done in the literature?

You should use a single temperature coupling group. I use the sd integrator (velocity Langevin dynamics) but I believe that you could also use velocity rescaling as your temperature coupling algorithm (I have not used it myself).

Usually with all bonds constrained by LINCS one uses a 2 fs timestep. I presume that you use a 0.5 fs timestep because that makes your conformation more stable? Again, this may just be hiding some problem with initial conformations and you may be able to move to a 2 fs timestep after some equilibration. As to the nstlist setting, it is common to update the neighbour list every 20 fs, but that will be dependent on your forcefield and your temperature. What have other people done with systems like this?

This will probably be my last comment on this thread. It is not usually possible to move from problems in system setup to a stable system with usable results within a week. I think that you would be best to aim for tracking down the problem and having a reliable simulation system by the time of your presentation.


-- original message --

Matthias Ernst Matthias.Ernst2 at student.kit.edu 
Thu Sep 6 17:55:43 CEST 2012
Previous message: [gmx-users] Effect of refcoord_scaling
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
> The problem could be gen_vel=no if you are not loading velocities that are consistent with your production temperature of 200 K. If, for example, you do not load any velocities at all, then the initial forces will quickly be scaled up to reach 200 K and this can cause large scale deformations/rearrangements.
I have generated velocities the step before (in NpT equilibration) and 
for the same temperature (200K). So, I would not think this causes the 
error as the should be transferred by using the resulting gro-file from 
NpT as input for production MD run.
> There are a number of other strange things (like nstlist=5 with dt=0.0005 -- nstlist seems unnecessarily frequent, the 200 K temperature, and your use of separate temperature coupling groups which does not give the correct ensemble). Still, those are tangents to your main question.
If I repeat it, I should do so with "correct" parameters, I think. Or do 
you suggest to repeat it with the very same parameter set as before?
To address the things you mentioned:
- What should or can nstlist be set to?
- 200K have been chosen because at higher temperatures, the test systems 
used to move even more.
- Concerning the coupling groups: you mean, I should use a group for 
protein and DNA as whole? But isn't restraining the DNA also a 
perturbation of the ensemble? At least, I thought that this way the 
temperature of the protein should not be affected by the restraints on 
the DNA. Tthis was the reason why I put several coupling groups and why 
I asked yesterday if it would be better not to couple the restrained DNA 
to the thermostate or set its temperature to 0.
> Nevertheless, if you don't have time to run any more simulations (as you stated earlier) then none of this is really going to help you. If, on the other hand, you do have time to run other simulations, then I still think you should start with repeating your result.
That's true. I fear, I will not be able to get it finished till next 
week but at least, I can try, start some jobs and see how far I can get.

Thank you very much for your help,

> Chris.
>   -- original message --
> I guess you want to see the production mdp file. Here it is:
> title                    = Production Simulation
> cpp                      = /lib/cpp
> define                   = -DPOSRES_DNA
> integrator               = md
> tinit                    = 0       ; Starting time
> dt                       = 0.0005 ; 2 femtosecond time step for integration
> nsteps                   = 100000000
> 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
> nstlist                  = 5
> ns-type                  = Grid
> pbc                      = xyz
> rlist                    = 0.9
> 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
> gen_vel                  = no
> constraints              = all-bonds
> constraint-algorithm     = Lincs
> unconstrained-start      = yes
> lincs-order              = 4
> lincs-iter               = 1
> lincs-warnangle          = 30
> Thanks for your help,
> Matthias

More information about the gromacs.org_gmx-users mailing list