Merge branch release-5-0 into release-5-1
[alexxy/gromacs.git] / src / gromacs / tools / compare.c
index 3b3947a6f19ff078bc375abfa3b10f4b7c16625b..5edc837bbc6c3f8ec7271e69380bf66b2ad5f344 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 /* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include <math.h>
 #include <stdio.h>
 #include <string.h>
-#include "main.h"
-#include "macros.h"
-#include "gromacs/utility/smalloc.h"
-#include "gromacs/fileio/futil.h"
-#include "sysstuff.h"
-#include "txtdump.h"
-#include "gmx_fatal.h"
-#include "names.h"
+
+#include "gromacs/fileio/enxio.h"
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/fileio/trxio.h"
-#include "gromacs/fileio/enxio.h"
-#include "mtop_util.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
 
 static void cmp_int(FILE *fp, const char *s, int index, int i1, int i2)
 {
@@ -89,11 +86,11 @@ static void cmp_us(FILE *fp, const char *s, int index, unsigned short i1, unsign
     {
         if (index != -1)
         {
-            fprintf(fp, "%s[%d] (%u - %u)\n", s, index, i1, i2);
+            fprintf(fp, "%s[%d] (%hu - %hu)\n", s, index, i1, i2);
         }
         else
         {
-            fprintf(fp, "%s (%u - %u)\n", s, i1, i2);
+            fprintf(fp, "%s (%hu - %hu)\n", s, i1, i2);
         }
     }
 }
@@ -478,6 +475,7 @@ static void cmp_top(FILE *fp, t_topology *t1, t_topology *t2, real ftol, real ab
         cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
         cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
         cmp_block(fp, &t1->mols, &t2->mols, "mols");
+        cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
         cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
     }
     else
@@ -524,6 +522,39 @@ static void cmp_groups(FILE *fp, gmx_groups_t *g0, gmx_groups_t *g1,
      */
 }
 
+static void cmp_rvecs_rmstol(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
+                             real ftol, real abstol)
+{
+    int    i, m;
+    double rms;
+
+    /* For a vector you are usally not interested in a relative difference
+     * on a component that is very small compared to the other components.
+     * Therefore we do the relative comparision relative to the RMS component.
+     */
+    rms = 0.0;
+    for (i = 0; (i < n); i++)
+    {
+        for (m = 0; m < DIM; m++)
+        {
+            rms += x1[i][m]*x1[i][m] + x2[i][m]*x2[i][m];
+        }
+    }
+    rms = sqrt(rms/(2*n*DIM));
+
+    /* Convert the relative tolerance into an absolute tolerance */
+    if (ftol*rms < abstol)
+    {
+        abstol = ftol*rms;
+    }
+
+    /* And now do the actual comparision */
+    for (i = 0; (i < n); i++)
+    {
+        cmp_rvec(fp, title, i, x1[i], x2[i], 0.0, abstol);
+    }
+}
+
 static void cmp_rvecs(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
                       gmx_bool bRMSD, real ftol, real abstol)
 {
@@ -545,54 +576,7 @@ static void cmp_rvecs(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
     }
     else
     {
-        for (i = 0; (i < n); i++)
-        {
-            cmp_rvec(fp, title, i, x1[i], x2[i], ftol, abstol);
-        }
-    }
-}
-
-
-/* Similar to cmp_rvecs, but this routine scales the allowed absolute tolerance
- * by the RMS of the force components of x1.
- */
-static void cmp_rvecs_rmstol(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
-                             real ftol, real abstol)
-{
-    int    i, m;
-    double d;
-    double ave_x1, rms_x1;
-
-    /* It is tricky to compare real values, in particular forces that
-     * are sums of lots of terms where the final value might be close to 0.0.
-     * To get a reference magnitude we calculate the RMS value of each
-     * component in x1, and then set the allowed absolute tolerance to the
-     * relative tolerance times this RMS magnitude.
-     */
-    ave_x1 = 0.0;
-    for (i = 0; i < n; i++)
-    {
-        for (m = 0; m < DIM; m++)
-        {
-            ave_x1 += x1[i][m];
-        }
-    }
-    ave_x1 /= n*DIM;
-
-    rms_x1 = 0.0;
-    for (i = 0; (i < n); i++)
-    {
-        for (m = 0; m < DIM; m++)
-        {
-            d       = x1[i][m] - ave_x1;
-            rms_x1 += d*d;
-        }
-    }
-    rms_x1 = sqrt(rms_x1/(DIM*n));
-    /* And now do the actual comparision with a hopefully realistic abstol. */
-    for (i = 0; (i < n); i++)
-    {
-        cmp_rvec(fp, title, i, x1[i], x2[i], ftol, abstol*rms_x1);
+        cmp_rvecs_rmstol(fp, title, n, x1, x2, ftol, abstol);
     }
 }
 
@@ -863,8 +847,8 @@ static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol,
     cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
     cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
 
-    cmp_int(fp, "inputrec->ePull", -1, ir1->ePull, ir2->ePull);
-    if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
+    cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
+    if (ir1->bPull && ir2->bPull)
     {
         cmp_pull(fp);
     }
@@ -915,7 +899,7 @@ static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol,
     cmp_cosines(fp, "et", ir1->et, ir2->et, ftol, abstol);
 }
 
-static void comp_pull_AB(FILE *fp, t_pull *pull, real ftol, real abstol)
+static void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol)
 {
     int i;
 
@@ -942,17 +926,17 @@ static void comp_state(t_state *st1, t_state *st2,
     if (st1->flags & (1<<estSVIR_PREV))
     {
         fprintf(stdout, "comparing shake vir_prev\n");
-        cmp_rvecs_rmstol(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, ftol, abstol);
+        cmp_rvecs(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, FALSE, ftol, abstol);
     }
     if (st1->flags & (1<<estFVIR_PREV))
     {
         fprintf(stdout, "comparing force vir_prev\n");
-        cmp_rvecs_rmstol(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, ftol, abstol);
+        cmp_rvecs(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, FALSE, ftol, abstol);
     }
     if (st1->flags & (1<<estPRES_PREV))
     {
         fprintf(stdout, "comparing prev_pres\n");
-        cmp_rvecs_rmstol(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, ftol, abstol);
+        cmp_rvecs(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, FALSE, ftol, abstol);
     }
     cmp_int(stdout, "ngtc", -1, st1->ngtc, st2->ngtc);
     cmp_int(stdout, "nhchainlength", -1, st1->nhchainlength, st2->nhchainlength);
@@ -1037,7 +1021,7 @@ void comp_tpx(const char *fn1, const char *fn2,
         }
         else
         {
-            if (ir[0].ePull != epullNO)
+            if (ir[0].bPull)
             {
                 comp_pull_AB(stdout, ir->pull, ftol, abstol);
             }
@@ -1094,14 +1078,7 @@ void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
     }
     if (cmp_bool(fp, "bF", -1, fr1->bF, fr2->bF))
     {
-        if (bRMSD)
-        {
-            cmp_rvecs(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, bRMSD, ftol, abstol);
-        }
-        else
-        {
-            cmp_rvecs_rmstol(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, ftol, abstol);
-        }
+        cmp_rvecs(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, bRMSD, ftol, abstol);
     }
     if (cmp_bool(fp, "bBox", -1, fr1->bBox, fr2->bBox))
     {