[gmx-users] Effect of refcoord_scaling

[Matthias Ernst]

> 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