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