# [gmx-users] is the tip4p MW velocity treated differently in gromacs 3 vs 4

Berk Hess gmx3 at hotmail.com
Wed Feb 4 11:19:25 CET 2009

```Hi,

First let me state again that these velocities are irrelevant for the dynamics,
because the velocities only enter the equation through momentum and kinetic energy
which have a factor mass, which is zero for vsites.

Looking at the code there are two issues.

One is with continuation (unconstrained start in 3.3).
In 3.3 unconstrained start would only remove the constraining of the initial
coordinates, but not the construction of vsites, which makes their velocity incorrect.
In 4.0 this is fully correct (but you need to start with correct x and v for the vsites).

Without continuation the velocities of vsites at step 0 and incorrect in both 3.3 and 4.0.
In 3.3 and 4.0 the velocities of vsites are determined from the differences in positions,
but the positions at step -1 for constraining the initial coordinates are not stored.
A clean solution (and more accurate) would be to determine the vsite velocities
directly with the vsite weights from the velocities of the masses, but that would
increase the computational cost.

A much simpler solution would be to simply always set all vsite velocities to zero.

I'll see if a can come up with a simple and clean solution.

Berk

From: gmx3 at hotmail.com
To: gmx-users at gromacs.org
Subject: RE: [gmx-users] is the tip4p MW velocity treated differently in	gromacs 3 vs 4
Date: Wed, 4 Feb 2009 09:03:29 +0100

Hi,

I always get confused by the atom names in TIP4P.
But MW is not a mass, but the vsite, right?
Velocities of vsites are completely irrelevant, except if you would want to analyze them.
I might have changed something to the initial velocities of vsites for 4.0.
I would guess that 4.0 has the correct initial vsite velocities,
but I would have to check to be sure.
Again, velocities of vsites do are completely irrelevant for anything,
except printing them in the output.

Berk

> Date: Wed, 4 Feb 2009 08:25:35 +0100
> From: spoel at xray.bmc.uu.se
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] is the tip4p MW velocity treated differently in	gromacs 3 vs 4
>
> Chris Neale wrote:
> > contrary to my previous stament, I am no longer convinced that it is
> > simply an initial constraints issue.
>
> There are two things to do:
> 1. run gmx3.3 and gmx4.0 with the 3.3 tpr file. This should give the
> same result in principle.
> 2. run gmx3.3 and gmx4.0 with their own tpr files. Here you can use
> gmxcheck -s1 -s2 to compare the files.
>
> The differences due to initial constraints are small usually, but they
> seem comparable in your case. Defining flexible should remove all
> difference due to constraints (unless constraints are turned on in the
> mdp file which selects LINCS over the default SETTLE).
>
> Anyway it is important to distinguish the apples from the oranges, hence
> point 1 is the best to test it. Obviously you are aware that the MW
> velocity does not in anyway influence the energy or the dynamics, since
> it merely shows the displacement of the particle divided by the time
> step. Therefore it would be good to compare energies as well from a 3.3
> run and a 4.0 run.
>
> >
> > In my last post, I stated: "when I define -DFLEXIBLE and make the tpr
> > using gromacs 3 and then run using gromacs 4, I do get gromacs 4 -
> > identical velocities.". However, thinking more on this has me believing
> > that this only indicates that the difference occurs at the mdrun stage.
> >
> > I did a few more tests to see what the difference was when turning on or
> > off the unconstrained start option
> > in either gromacs 3 or gromacs 4. I see a very small diffence in
> > velocities compared to the MW difference between
> > gromacs 3 or gromacs 4 mdruns.
> >
> > Below showing gromacs 4, using unconstrained_start=no vs.
> > unconstrained_start=yes has differences, but very minor.
> >
> > diff gmx4.0.3_noUCS/feoff.gro gmx4.0.3/feoff.gro ...
> > <    51SOL     OW 1151   6.816   4.071   1.558  0.0549  0.4777 -0.4315
> > <    51SOL    HW1 1152   6.851   4.082   1.470 -0.1304  0.6992 -0.4795
> > <    51SOL    HW2 1153   6.727   4.106   1.553 -0.2840 -0.3333 -0.3852
> > <    52SOL     OW 1155   3.564   2.550   1.079  0.0951  0.0648  0.1748
> > <    52SOL    HW1 1156   3.521   2.598   1.009  1.5023  0.4847 -0.4757
> > <    52SOL    HW2 1157   3.551   2.604   1.157  0.0156  0.5140 -0.1424
> > ---
> >>    51SOL     OW 1151   6.816   4.071   1.558  0.0554  0.4775 -0.4315
> >>    51SOL    HW1 1152   6.851   4.082   1.470 -0.1306  0.6992 -0.4793
> >>    51SOL    HW2 1153   6.727   4.106   1.553 -0.2903 -0.3308 -0.3857
> >>    52SOL     OW 1155   3.564   2.550   1.079  0.0951  0.0648  0.1747
> >>    52SOL    HW1 1156   3.521   2.598   1.009  1.5034  0.4843 -0.4719
> >>    52SOL    HW2 1157   3.551   2.604   1.157  0.0147  0.5151 -0.1438
> > ...
> >
> > with some, but not all MW showing small differences:
> > <    54SOL     OW 1163   5.835   7.123   3.589 -0.0214 -0.1603  0.5173
> > <    54SOL    HW1 1164   5.788   7.068   3.651  0.4698 -0.5196  0.5833
> > <    54SOL    HW2 1165   5.841   7.208   3.632 -0.8849 -0.0052  0.3674
> > <    54SOL     MW 1166   5.830   7.127   3.602 -0.0578 -0.2071  0.4856
> > ---
> >>    54SOL     OW 1163   5.835   7.123   3.589 -0.0214 -0.1604  0.5172
> >>    54SOL    HW1 1164   5.788   7.068   3.651  0.4694 -0.5205  0.5835
> >>    54SOL    HW2 1165   5.841   7.208   3.632 -0.8844 -0.0018  0.3686
> >>    54SOL     MW 1166   5.829   7.127   3.602 -0.0578 -0.2071  0.4856
> >
> >
> > Below showing gromacs 3, using unconstrained_start=no vs.
> > unconstrained_start=yes has differences, but very minor.
> >
> > diff gmx3.3.1_noUCS/feoff.gro gmx3.3.1/feoff.gro ...
> > <    51SOL     OW 1151   6.816   4.071   1.558  0.0550  0.4778 -0.4315
> > <    51SOL    HW1 1152   6.851   4.082   1.470 -0.1304  0.6992 -0.4795
> > <    51SOL    HW2 1153   6.727   4.106   1.553 -0.2841 -0.3333 -0.3852
> > <    51SOL     MW 1154   6.809   4.077   1.546  0.0254 -0.0273  0.0172
> > <    52SOL     OW 1155   3.564   2.550   1.079  0.0951  0.0648  0.1749
> > <    52SOL    HW1 1156   3.521   2.598   1.009  1.5023  0.4846 -0.4757
> > <    52SOL    HW2 1157   3.551   2.604   1.157  0.0157  0.5140 -0.1425
> > <    52SOL     MW 1158   3.557   2.563   1.080 -0.0486  0.0225 -0.0022
> > ---
> >>    51SOL     OW 1151   6.816   4.071   1.558  0.0554  0.4775 -0.4315
> >>    51SOL    HW1 1152   6.851   4.082   1.470 -0.1306  0.6992 -0.4793
> >>    51SOL    HW2 1153   6.727   4.106   1.553 -0.2903 -0.3308 -0.3857
> >>    51SOL     MW 1154   6.809   4.077   1.546  0.0219 -0.0279  0.0237
> >>    52SOL     OW 1155   3.564   2.550   1.079  0.0951  0.0648  0.1747
> >>    52SOL    HW1 1156   3.521   2.598   1.009  1.5034  0.4843 -0.4719
> >>    52SOL    HW2 1157   3.551   2.604   1.157  0.0147  0.5151 -0.1438
> >>    52SOL     MW 1158   3.557   2.563   1.080 -0.0421  0.0144  0.0060
> >
> > ...
> >
> > <    54SOL     OW 1163   5.835   7.123   3.589 -0.0215 -0.1603  0.5173
> > <    54SOL    HW1 1164   5.788   7.068   3.651  0.4698 -0.5196  0.5834
> > <    54SOL    HW2 1165   5.841   7.208   3.632 -0.8849 -0.0051  0.3674
> > <    54SOL     MW 1166   5.830   7.127   3.602  0.1837 -0.0398  0.1193
> > ---
> >>    54SOL     OW 1163   5.835   7.123   3.589 -0.0214 -0.1604  0.5172
> >>    54SOL    HW1 1164   5.788   7.068   3.651  0.4694 -0.5205  0.5835
> >>    54SOL    HW2 1165   5.841   7.208   3.632 -0.8844 -0.0018  0.3686
> >>    54SOL     MW 1166   5.830   7.127   3.602  0.1879 -0.0399  0.1103
> > ...
> >
> > In the above cases, other specified .mdp options were:
> >
> > integrator          =  md
> > comm_mode           =  linear
> > nstcomm             =  1
> > comm_grps           =  System
> > nstlist             =  5
> > ns_type             =  grid
> > pbc                 =  xyz
> > coulombtype         =  PME
> > rcoulomb            =  0.9
> > fourierspacing      =  0.12
> > pme_order           =  4
> > vdwtype             =  cut-off
> > rvdw_switch         =  0
> > rvdw                =  1.4
> > rlist               =  0.9
> > constraints         =  none
> > energygrps          =  System
> > nsteps              =  0
> > dt                  =  0.004
> > gen_vel             =  no
> > unconstrained-start =  <yes or no depending on the test>
> >
> > ######
> >
> > And note that if I set
> > constraints         =  all-bonds
> > constraint_algorithm=  lincs
> > lincs-iter          =  1
> > lincs-order         =  6
> >
> > then the difference between unconstrained-start=yes or =no was similarly
> > small.
> >
> > Therefore it seems to me that something other than initial constraints
> > is leading to this difference.
> >
> > Chris.
> >
> > -- original message --
> >
> > Thanks David, you are correct about initial step constraints. More
> > information inline below, but my question has been resolved.
> >
> >  >Chris Neale wrote:
> >  >> Hello,
> >  >>
> >  >> Does anybody know if there is a reason why the .gro output velocities
> >  >> would be different for tip4p MW in a zero-step mdrun between gromacs 3
> >  >> and gromacs 4 (3.3.1 and 3.3.3 are the same, and are different from
> >  >> 4.0.2 and 4.0.3, which are themselves the same).
> >  >Is this with the same tpr?
> >
> > This was utilizing a tpr generated by the same version of gromacs as was
> > used from the mdrun (3.3.1 grompp for 3.3.1 mdrun, etc.) although the
> > input was identical (with the exception of any undefined .mdp options
> > whose default values changes between v3 and v4, although I can't think
> > of one that affects only MW).
> >
> > Based on your question, I did redo the 3.3.1 grompp to generate a single
> > .tpr that was used for both 3.3.1 mdrun and 4.0.3 mdrun. The results are
> > the same:
> >
> > \$ diff 4.gro 3.gro | head
> > 1156c1156
> > <    51SOL     MW 1154   6.809   4.077   1.546  0.0120  0.3740 -0.4572
> > ---
> >  >    51SOL     MW 1154   6.809   4.077   1.546  0.0219 -0.0279  0.0237
> > 1160c1160
> > <    52SOL     MW 1158   3.557   2.563   1.080  0.2917  0.2007  0.0277
> > ---
> >  >    52SOL     MW 1158   3.557   2.563   1.080 -0.0421  0.0144  0.0060
> >
> >  >
> >  >It could have to do with initial step constraints.
> >
> > I had considered that, but my .mdp uses unconstrained-start =  yes.
> >
> > Based on your suggestion I have attempted setting:
> > 1. constraints = none (get the same difference)
> > 2. also define = -DFLEXIBLE (get the same difference)
> >
> > However, when I define -DFLEXIBLE and make the tpr using gromacs 3 and
> > then run using gromacs 4
> > I do get gromacs 4 - identical velocities.
> >
> > \$ diff 4from3.gro 4.gro | head
> > \$ diff 4from3.gro 3.gro | head
> > 1156c1156
> > <    51SOL     MW 1154   6.809   4.077   1.546  0.0120  0.3740 -0.4572
> > ---
> >  >    51SOL     MW 1154   6.809   4.077   1.546  0.0219 -0.0279  0.0237
> > 1160c1160
> > <    52SOL     MW 1158   3.557   2.563   1.080  0.2917  0.2007  0.0277
> > ---
> >  >    52SOL     MW 1158   3.557   2.563   1.080 -0.0421  0.0144  0.0060
> > 1164c1164
> > <    53SOL     MW 1162   2.458   4.404   3.060 -0.0763  0.0165 -0.5712
> >
> > So thank you, this explains it.
> >
> > Chris.
> >
> >  >>
> >  >> diff gmx4.0.3/feoff.gro gmx3.3.1/feoff.gro | head
> >  >> 1156c1156
> >  >> <    51SOL     MW 1154   6.809   4.077   1.546  0.0120  0.3740 -0.4572
> >  >> ---
> >  >>  >    51SOL     MW 1154   6.809   4.077   1.546  0.0219 -0.0279  0.0237
> >  >> 1160c1160
> >  >> <    52SOL     MW 1158   3.557   2.563   1.080  0.2917  0.2007  0.0277
> >  >> ---
> >  >>  >    52SOL     MW 1158   3.557   2.563   1.080 -0.0421  0.0144  0.0060
> >  >> ...
> >  >> (continues for all MW)
> >  >>
> >  >> This difference does not appear to be related to solvent optimization:
> >  >>
> >  >> \$ diff gmx4.0.3/feoff.gro gmx4.0.3_GMX_NO_SOLV_OPT/feoff.gro
> >  >>
> >  >> I might just chalk this up to different mdp option defaults, but it
> >  >> seems strange that there is a difference for all MW, but not for any
> >  >> other atoms in the system.
> >  >>
> >  >> \$ cat dpc50_md11.mdp
> >  >> integrator          =  md
> >  >> comm_mode           =  linear
> >  >> nstcomm             =  1
> >  >> comm_grps           =  System
> >  >> nstlist             =  5
> >  >> ns_type             =  grid
> >  >> pbc                 =  xyz
> >  >> coulombtype         =  PME
> >  >> rcoulomb            =  0.9
> >  >> fourierspacing      =  0.12
> >  >> pme_order           =  4
> >  >> vdwtype             =  cut-off
> >  >> rvdw_switch         =  0
> >  >> rvdw                =  1.4
> >  >> rlist               =  0.9
> >  >> constraints         =  all-bonds
> >  >> constraint_algorithm=  lincs
> >  >> lincs-iter          =  1
> >  >> lincs-order         =  6
> >  >> energygrps          =  System
> >  >> nsteps              =  0
> >  >> dt                  =  0.004
> >  >> gen_vel             =  no
> >  >> unconstrained-start =  yes
> >  >>
> >  >> Note that I explicitly included the tip4p.itp from gromacs 3.3.1 in
> > both
> >  >> runs.
> >  >>
> >  >> Perhaps this is related to the gromacs 3 vsite problem?
> >  >> http://www.gromacs.org/pipermail/gmx-users/2008-November/037659.html
> >  >>
> >  >> Thanks,
> >  >> Chris.
> >
> > _______________________________________________
> > gmx-users mailing list    gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at http://www.gromacs.org/search before posting!
> > Please don't post (un)subscribe requests to the list. Use the www
> > interface or send it to gmx-users-request at gromacs.org.
> > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
> spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.