[gmx-users] shake for water

[David van der Spoel]

Berk Hess wrote:
> 
> 
>  > Date: Thu, 5 Feb 2009 19:35:09 +0100
>  > From: spoel at xray.bmc.uu.se
>  > To: gmx-users at gromacs.org
>  > Subject: Re: [gmx-users] shake for water
>  >
>  > David Mobley wrote:
>  > > All,
>  > >
>  > > A quick question on constraints... I'm using TIP4P-Ew in gromacs 3.3.2
>  > > and am concerned with reproducing energies from another code very
>  > > precisely for several specific snapshots. I am doing a zero-step mdrun
>  > > of a setup with one small molecule and two tip4p-ew water molecules.
>  > >
>  > > Anyway, I have set the shake tolerance to 1e-12 in the mdp file, but
>  > > to my surprise the internal water distances are good only to 1e-06 and
>  > > 1e-07. Is this expected behavior? Note that I am running in double
>  > > precision. I assumed that, er, the distances should converge to the
>  > > shake tolerance.
>  > Well, the documentation might be lacking, but the code tells the truth.
>  > It seems that the tolerance is used on the distance squared, which is
>  > consistent with your observation of a precision of 1e-6. So try 1e-24.
> 
> No, it is not the square.
> The code does a small 'a' approximation: (1+a)^2=1+2a+a^2 is approx 1+2a.
> I have also tested/benchmark shake in gromacs for my first lincs paper
> and the plincs paper and it always behaved the way I thought it would.


First we compute the inverse square of the shake distance dA in tt[ll] 
for each shake pair:

     if (bFEP)
       toler = sqr(L1*ip[type].shake.dA + lambda*ip[type].shake.dB);
     else
       toler = sqr(ip[type].shake.dA);
     dist2[ll] = toler;
     tt[ll] = 1.0/(toler*tol2);
   }

Then in the shake iteration we compute the difference between the 
squared distances (variable diff below):

       tx      = xp[ix]-xp[jx];
       ty      = xp[iy]-xp[jy];
       tz      = xp[iz]-xp[jz];
       rpij2   = tx*tx+ty*ty+tz*tz;
       toler   = dist2[ll];
       diff    = toler-rpij2;

Now we multiply the diff with tt[ll], in other words we get
iconv = (1-rpij2/dist2)/(2 tol)

       /* iconv is zero when the error is smaller than a bound */
       iconv   = fabs(diff)*tt[ll];

In other words, the tolerance operates on the squared distance.




> 
> The problem is not that you need more iterations than 1000?
> 
> And why not use settle that does it right at full precision at once?
> 
> Berk
> 
> 
>  >
>  > >
>  > > Thanks,
>  > > David
>  > > _______________________________________________
>  > > 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.
>  > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> 
> ------------------------------------------------------------------------
> Express yourself instantly with MSN Messenger! MSN Messenger 
> <http://clk.atdmt.com/AVE/go/onm00200471ave/direct/01/>
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> 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