[gmx-users] Conserved Energy in csvr NVT runs

Thu Apr 17 17:04:09 CEST 2014

Dear Mark, Dear Users,

Thank you for the detailed and prompt answer!
I'm sorry I probably wasn't clear enough, and I didn't attach the full
input.  So I'll attach to avoid confusion.

>IIRC there is a conserved quantity written, but as you can see in Bussi's
>original paper, it has a t dependence. This has been discussed on this list
>before, if you search. tau-t probably affects the slope.

If you refer to the J. Chem. Phys. 126, 014101 (2007) paper, you will
exactly see there, what I wish to see in my simulation.
On Fig 2, they plotted the Conserved Quantity ($\widetilede{H}_{NVT}$) with
respect time for  a LJ system.
With 5 fs timestep they have nice flat, non-drifting oscillation.  I would
like to check the same quantity to set my timestep for a proper value on
the timescale I do my simultaion.

In my simulation, I printed  the "Conserved-En." by g_energy_d for the
timesteps:  0.5 fs , 0.25 and  0.1 fs when I applied  md-vv integrator (As
I have h-bond constraint, I think 1.0 fs should be OK, 0.5 -0.25 fs safe,
0.1 fs too careful ).
With comm-mode=none, integrator=md-vv I get in the two first case clear
drift and in case of 0.1 fs I get proper conservation.
I'm happy with that, but with constrained h-bond length I'm a bit
surprised, anyway in NVE I didn't need such small time step.

I also checked integrator=md, and it gives much smaller drift for any
timestep. See the attached drift graph for different timesteps and
integrators.

I got the hint from a friend, that I should use sd instead of md-vv or md
integrator.
If I do so, "Conserved-En" disappears from the options in g_energy_d. Total
energy remains nicely oscillating.

So the question now is: What is the proper integrator (and other settings)
to use, and which output energy type (if any) equals to
$\widetilede{H}_{NVT}$, defined by eq 15 in Giovanni's paper?

If the  $\widetilede{H}_{NVT}$ is printed as "Conserved-En", and md is the
proper integrator, is this small timestep (0.25-0.1 fs) a normal value for
bond constrained system?

>You haven't said, but I assume you have a dipeptide in vacuum, which I
>think should need angular if the action of the thermostat can add kinetic
>energy to the COM dof.

Yes, and sorry again to keep this information. I try to make non-periodic
simulation in vacuum.
To understand perfectly could you tell me how many degrees of freedom left
with  pbc=no option?
a) 3*N-Nconst-0  # I would use: comm-mode = none
b) 3*N-Nconst-3  # translation is not thermostated; I would use: comm-mode
= linear
c) 3*N-Nconst-6  # translation and rotation is not thermostated; I would
use: comm-mode = angular

(N=number of atoms, Nconst=number of constraints)
I have seen in other codes the c) option for initial velocity generation,
what is the case with gmx? What is the right option for comm-mode?

Any suggestion and comment concerting my input is greatly acknowledged!

Thank you in advance,

Janos

2014-04-16 21:18 GMT+02:00 Mark Abraham <mark.j.abraham at gmail.com>:

> On Apr 16, 2014 2:04 PM, "janos.daru" <janos.daru at gmail.com> wrote:
> >
> > Dear All,
> >
> > I run simple md job (GMX version:4.6.5 double precision) for
> > alanine-dipeptid with careful options.
> > nstlist                  = 1
> > ns_type               = simple
> > pbc                     = no
> > rlist                     = 0
> > epsilon_r                = 1
> > coulombtype          = cut-off
> > rcoulomb                 = 0
> > vdwtype                  = cut-off
> > rvdw                     = 0
> > dispcorr                 = no
> > constraints              = h-bonds
> > constraint-algorithm     = Lincs
> > lincs-order              = 4
> > lincs-iter               = 2
> >
> > For NVE runs I get very nice conserved energy values: ~10E-5
> > K/numberOfDegreesOfFreedom/ps
> >
> > But when I change to v-rescale with the following settings, I get clear
> > drift(~1 K/numberOfDegreesOfFreedom/ps).
> >
> > tcoupl                   = v-rescale
> > tc-grps                  = Protein
> > tau_t                    = 0.05
> > ref_t                    = 300
> > nsttcouple               = 1
> >
> > As total energy and temperature makes a proper oscillation I have the
> > filling, that "Conserved Energy" is not calculated as I expected, or I
> miss
> > something important form the input.
>
> IIRC there is a conserved quantity written, but as you can see in Bussi's
> original paper, it has a t dependence. This has been discussed on this list
> before, if you search. tau-t probably affects the slope.
>
> > An other observation is that in NVE runs I need to use comm_mode=none to
> get
> > nice energy conservation, while in NVT I get slightly better results with
> > comm_mode=angular for longer runs (aprox 2.5 ns).
>
> You haven't said, but I assume you have a dipeptide in vacuum, which I
> think should need angular if the action of the thermostat can add kinetic
> energy to the COM dof.
>
> > Could some experienced user comment on these observations and give me
> > advice, how to proceed to get proper NVT ensamble without artifacts of
> the
> > drift?
>
> You probably did. What is your time step?
>
> > I didn't find in the manual options for global and local thermostats. Can
> I
> > use local (each degrees of freedom) thermostat in gromacs?
>
> Not really. Best you can do is a temperature coupling group per atom.
>
> Mark
>
> > Thank you in advance,
> >
> > Janos
> >
> >
> >
> >
> >
> > --
> > View this message in context:
>
> http://gromacs.5086.x6.nabble.com/Conserved-Energy-in-csvr-NVT-runs-tp5015881.html
> > Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>