From 6a602ac85f08cf82d138038d9960001e4040ffff Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Wed, 24 Aug 2016 11:00:19 +0200
Subject: [PATCH] Reduce rounding errors in SETTLE

The parameters for SETTLE are now computed in double precision.
This lowers the systematic error in settle.

Change-Id: I537a618830ab149d5f33251d9edca3aa9c43404c
---
 src/gromacs/mdlib/csettle.cpp | 55 ++++++++++++++++++++++++-------------------
 1 file changed, 31 insertions(+), 24 deletions(-)

diff --git a/src/gromacs/mdlib/csettle.cpp b/src/gromacs/mdlib/csettle.cpp
index 092ee94ac9..082b14e176 100644
--- a/src/gromacs/mdlib/csettle.cpp
+++ b/src/gromacs/mdlib/csettle.cpp
@@ -94,44 +94,43 @@ typedef struct gmx_settledata
 } t_gmx_settledata;
 
 
-static void init_proj_matrix(settleparam_t *p,
-                             real invmO, real invmH, real dOH, real dHH)
+static void init_proj_matrix(real invmO, real invmH, real dOH, real dHH,
+                             matrix inverseCouplingMatrix)
 {
-    real   imOn, imHn;
-    matrix mat;
-
-    p->imO = invmO;
-    p->imH = invmH;
-    /* We normalize the inverse masses with imO for the matrix inversion.
+    /* We normalize the inverse masses with invmO for the matrix inversion.
      * so we can keep using masses of almost zero for frozen particles,
      * without running out of the float range in invertMatrix.
      */
-    imOn = 1;
-    imHn = p->imH/p->imO;
+    double invmORelative = 1.0;
+    double invmHRelative = invmH/static_cast<double>(invmO);
+    double distanceRatio = dHH/static_cast<double>(dOH);
 
     /* Construct the constraint coupling matrix */
-    mat[0][0] = imOn + imHn;
-    mat[0][1] = imOn*(1 - 0.5*dHH*dHH/(dOH*dOH));
-    mat[0][2] = imHn*0.5*dHH/dOH;
+    matrix mat;
+    mat[0][0] = invmORelative + invmHRelative;
+    mat[0][1] = invmORelative*(1.0 - 0.5*gmx::square(distanceRatio));
+    mat[0][2] = invmHRelative*0.5*distanceRatio;
     mat[1][1] = mat[0][0];
     mat[1][2] = mat[0][2];
-    mat[2][2] = imHn + imHn;
+    mat[2][2] = invmHRelative + invmHRelative;
     mat[1][0] = mat[0][1];
     mat[2][0] = mat[0][2];
     mat[2][1] = mat[1][2];
 
-    gmx::invertMatrix(mat, p->invmat);
-
-    msmul(p->invmat, 1/p->imO, p->invmat);
+    invertMatrix(mat, inverseCouplingMatrix);
 
-    p->invdOH = 1/dOH;
-    p->invdHH = 1/dHH;
+    msmul(inverseCouplingMatrix, 1/invmO, inverseCouplingMatrix);
 }
 
 static void settleparam_init(settleparam_t *p,
                              real mO, real mH, real invmO, real invmH,
                              real dOH, real dHH)
 {
+    /* We calculate parameters in double precision to minimize errors.
+     * The velocity correction applied during SETTLE coordinate constraining
+     * introduces a systematic error of approximately 1 bit per atom,
+     * depending on what the compiler does with the code.
+     */
     double wohh;
 
     p->mO     = mO;
@@ -140,13 +139,21 @@ static void settleparam_init(settleparam_t *p,
     p->wh     = mH/wohh;
     p->dOH    = dOH;
     p->dHH    = dHH;
-    p->rc     = dHH/2.0;
-    p->ra     = 2.0*mH*sqrt(dOH*dOH - p->rc*p->rc)/wohh;
-    p->rb     = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
+    double rc = dHH/2.0;
+    double ra = 2.0*mH*std::sqrt(dOH*dOH - rc*rc)/wohh;
+    p->rb     = std::sqrt(dOH*dOH - rc*rc) - ra;
+    p->rc     = rc;
+    p->ra     = ra;
     p->irc2   = 1.0/dHH;
 
-    /* For projection: connection matrix inversion */
-    init_proj_matrix(p, invmO, invmH, dOH, dHH);
+    /* For projection: inverse masses and coupling matrix inversion */
+    p->imO    = invmO;
+    p->imH    = invmH;
+
+    p->invdOH = 1.0/dOH;
+    p->invdHH = 1.0/dHH;
+
+    init_proj_matrix(invmO, invmH, dOH, dHH, p->invmat);
 
     if (debug)
     {
-- 
2.11.4.GIT