[gmx-users] program crash

[Mark Abraham]

> Dear GMX users,
>
> My program crashed very soon after it started. I had done EM, and the
> conformation looked fine in VMD. When I investigated what's wrong, I
> outputted the coordinates, forces, velocities and log file in every
> step. I set some parameters in mdp file as follows:
>
> gen_vel                  = no; yes
> gen_temp                 = 100
> gen_seed                 = 1223
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling
> Tcoupl                   = berendsen  ; nose-hoover ;berendsen
> ; Groups to couple separately
> tc-grps                  = protein
> ; Time constant (ps) and reference temperature (K)
> tau_t                    = 0.1
> ref_t                    = 100
> ; Pressure coupling
> Pcoupl                   =  no;
>
>
> The following is from my log file:
>
> There are 1802 atoms in your xtc output selection
>
> ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
> H. J. C. Berendsen, J. P. M. Postma, A. DiNola and J. R. Haak
> Molecular dynamics with coupling to an external bath
> J. Chem. Phys. 81 (1984) pp. 3684-3690
> -------- -------- --- Thank You --- -------- --------
>
>    Energies (kJ/mol)
>            Bond          Angle       G96Angle Ryckaert-Bell.
> LJ-14
>     2.23601e+03    3.10376e+02    9.73565e+02    4.90946e+03    3.92667e
> +04
>      Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic
> En.
>     3.73370e+03   -3.85288e+03   -3.22634e+03    4.43506e+04    1.21982e
> +04
>    Total Energy    Temperature Pressure (bar)
>     5.65488e+04    5.43368e+02    0.00000e+00
>
>            Step           Time         Lambda
>               1        0.00200        0.00000
>
>
> This is the first run after EM. I set "gen_vel = no" on purpose. When I
> checked the trr file, the initial velocities were zeros, so that's
> right. However, the log file showed that the initial KE and T were not
> zeros. This confuses me. How does gromacs calculate KE and T? Another
> strange thing is that the initial temperature is 500K, not 100K even
> though I set ref_t =100. Any help will be highly appreciated!

As you'll know from undergrad courses, you can't have temperature without
kinetic energy. So if you generate no velocities, your initial temperature
need have no relationship to the one you "set" with ref_t - remember that
that is a parameter for a negative-feedback algorithm and takes time to
achieve equilibrium. I suggest you read the relevant manual sections
closely. I can't explain where the KE is coming from, unless there's some
force integrations "before" step 1 such that atoms have started moving.

In any case, there is no reason to want gen_vel = no unless you are
keeping velocities from a previous simulation.

Mark