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

David van der Spoel spoel at xray.bmc.uu.se
Wed Feb 4 08:25:35 CET 2009


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



More information about the gromacs.org_gmx-users mailing list