[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.
Chris.
-- 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,
Matthias
> Chris.
>
> -- original message --
>
> 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
More information about the gromacs.org_gmx-users
mailing list