[gmx-users] System does not heat up to 300K.

[Justin Lemkul]

On 6/11/13 2:07 PM, Marc Hömberger wrote:
> Hi,
>
> I have a system setup (minimized etc.) and wanted to heat this system
> gradually with simulated annealing. I set up the mdp file as you can see
> below (simulated annealing part is marked bold). But when the heating is
> finished, the log file tells me that the temperature of the system is only
> 1K (see below).
>
> Can someone tell me what I am doing wrong?
>

You've specified heating to 300K over 300,000 ps (300 ns), thus 1 K per 1 ns. 
Your simulation calls for 1 ns of time, thus 1 K should be achieved.

-Justin

>
> Best,
> Marc
>
> ;************ PREPROCESSING *********************************************
> ;include                =
> ;define                 =
> ;************ RUN CONTROL ***********************************************
> dt                      = 0.002         ; timestep, 2 fs
> integrator              = md            ; leap-frog integrator
> nsteps                  = 500000     ; 2 * 500000 = 1000 ps
> tinit                   = 0             ; starting time, ps
> init_step               = 0             ; starting step
> comm-mode               = linear        ; remove COM translation
> nstcomm                 = 10            ; freq steps comm-mode
> comm-grps               = system        ; groups COM removal ; protein
> nstcalcenergy           = -1            ; freq calc ener; default = nstlist
> ;************ DYNAMICS **************************************************
> ;### Initial velocity assignment
> gen_vel                 = yes           ; assign velocities from Maxwell
> distribution
> gen_temp                = 300           ; temperature for Maxwell
> distribution
> gen_seed                = -1            ; random generator seed; -1 from
> job ID
> ;************ OUTPUT CONTROL ********************************************
> ;### xtc
> nstxtcout               = 2500          ; .xtc freq write coor
> xtc-precision           = 1000          ; .xtc precision
> xtc_grps                = system        ; .xtc saved groups
> ;### trr
> nstxout                 = 500000        ; .trr freq write coor
> nstvout                 = 500000        ; .trr freq write vel
> nstfout                 = 0             ; .trr freq write forces
> ;### log
> nstlog                  = 2500          ; .log freq write ener
> ;### edr
> nstenergy               = 2500          ; .edr freq write ener; =
> i*nstcalcenergy
> energygrps              = protein non-protein   ; .edr saved groups
> ;************ NON-BONDED INTERACTIONS ***********************************
> ns_type                 = grid          ; search neighbor grid cells
> nstlist                 = 10            ; freq update NB list
>
> vdwtype                 = switch        ; switching function for VdW
> coulombtype             = PME           ; PME for long-range electrostatics
> rvdw_switch             = 0.9           ; switching dist VdW
> rvdw                    = 1.0           ; cutoff VdW
> rcoulomb                = 1.2           ; cutoff elec
> rlist                   = 1.2           ; cutoff NB list; must = rcoulomb
> when PME
>
> epsilon-r               = 1             ; dielectric const
> DispCorr                = EnerPres      ; Long-range VdW dispersion
> correction
>
> ;### PME & Ewald
> ewald_rtol              = 1.0e-5        ; direct space tolerance;
> 1.0e-6/2.2e-5
> fourierspacing          = 0.12          ; max grid spacing for FFT
> pme_order               = 6             ; interpolation order; depends on
> cutoff
> ewald_geometry          = 3d            ; Ewald in 3d
> epsilon_surface         = 0             ; dipole correction for 3d Ewald
> optimize_fft            = yes           ; start-up FFT grid optimization
> ;************ BOUNDARY CONTROL ******************************************
> pbc                     = xyz           ; 3-D PBC
> ;************ CONSTRAINTS CONTROL ***************************************
> constraints             = all-bonds     ; all bonds (even heavy atom-H
> bonds) constrained
> continuation            = no            ; constraints applied to initial
> conf
> constraint_algorithm = lincs    ; holonomic constraints
> lincs_iter              = 1             ; accuracy of LINCS
> lincs_order             = 4             ; also related to accuracy
> lincs_warnangle         = 30            ; max allowed bond angle before
> warning
> ;************ TEMPERATURE CONTROL ***************************************
> tcoupl                  = nose-hoover   ; Nose-Hoover T-coupling; v-rescale
> tc-grps                 = system        ; two coupling groups - more
> accurate
> tau_t                   = 0.6           ; coupling time constant at equil,
> in ps
> ref_t                   = 300           ; reference temperature, one for
> each group, in K
> ;************ PRESSURE CONTROL ******************************************
> pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in
> NPT
> pcoupltype              = isotropic     ; uniform scaling of box vectors
> tau_p                   = 1.0           ; coupling time constant at equil,
> in ps
> ref_p                   = 1.0           ; reference pressure, in bar
> compressibility         = 4.5e-5        ; isothermal compressibility of
> water at 1 atm, bar^-1
> ;************ SIMULATED ANNEALING ******************************************
> *annealing               = single        ; type of annealing (single
> sequence of annealing points)*
> *annealing_npoints       = 2             ; # of temperature points*
> *annealing_time          = 0 300000      ; time window in which temperature
> shall be increased*
> *annealing_temp          = 0 300         ; temperature at timepoints*
>
>
> LOG-FILE:
> Writing checkpoint, step 500000 at Thu Jun  6 21:23:55 2013
>
>
> *Current ref_t for group System:      1.0*
>     Energies (kJ/mol)
>              U-B    Proper Dih.  Improper Dih.      CMAP Dih.          LJ-14
>      8.64012e+03    5.98740e+03    4.75546e+02   -1.72734e+03    4.25420e+03
>       Coulomb-14        LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.
>      4.78689e+04    1.75858e+05   -5.83865e+03   -1.21361e+06   -1.12292e+05
>   Position Rest.      Potential    Kinetic En.   Total Energy    Temperature
>      7.52594e+00   -1.09037e+06    1.35115e+05   -9.55257e+05    2.17729e+02
>   Pres. DC (bar) Pressure (bar)   Constr. rmsd
>     -1.35977e+02    1.38325e+02    2.75899e-05
>

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================