Merge branch release-2019

Change-Id: I71b540586e32f54f3536c5a59eb99d89ad0b3639

diff --git a/docs/release-notes/2019/2019.4.rst b/docs/release-notes/2019/2019.4.rst
index caa2a0ab6403e2478f3950f7e7d215da0b37397f..3d3f59af7f316d73185c903404cccf40f2214c6d 100644 (file)
--- a/docs/release-notes/2019/2019.4.rst
+++ b/docs/release-notes/2019/2019.4.rst
@@ -16,6 +16,18 @@ in the :ref:`release-notes`.
 Fixes where mdrun could behave incorrectly
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+Fix (unlikely) missing bonded forces with CUDA GPUs and domain decomposition
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Forces could be missing for bonded interactions computed on CUDA GPUs with
+domain decomposition when there are non-local bonded interactions, but no
+non-local non-bonded interactions between two domains. Note that this is
+extremely unlikely to happen, since the distance between the bonded atoms
+needs to be larger than the pair-list cut-off distance and there should be no
+other non-local atoms within the pair-list cut-off distance.
+
+:issue:`3063`
+
 Fix segmentation fault in grompp and mdrun with cosine COM pulling
 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
@@ -32,7 +44,6 @@ The tool would fail when not being provided with a library file to read in.
 
 :issue:`3012`
 
-
 Fix bug in gmx anaeig
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
@@ -42,6 +53,18 @@ of eigenvectors was less than the three times the number of atoms.
 
 :issue:`2972`
 
+Fix issue with demux.pl script
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The trajectories could become discontinuous with simulations longer than 100ns
+and exchange strides that are not a multiple of 1 ps. This only affected the
+post-processing of trajectories generated from replica exchange simulations.
+
+Made gmx disre work with non-consecutively labeled restraints
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2953`
+
 Fixes that affect portability
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/scripts/demux.pl b/scripts/demux.pl
index f8e2f558ff8c63b1b4ba3cb47ad6ef501dd9bdda..23cb6ef52bdeee392f6fb58959a369b95f0a725a 100755 (executable)
--- a/scripts/demux.pl
+++ b/scripts/demux.pl
@@ -30,7 +30,7 @@ open (TEMP,">$temp") || die("Opening $temp for writing");
 sub pr_order {
     my $t     = shift;
     my $nrepl = shift;
-    printf(NDX "%-20g",$t);
+    printf(NDX "%-20.2f",$t);
     for(my $k=0; ($k<$nrepl); $k++) {
        my $oo = shift;
        printf(NDX "  %3d",$oo);
@@ -41,7 +41,7 @@ sub pr_order {
 sub pr_revorder {
     my $t     = shift;
     my $nrepl = shift;
-    printf(TEMP "%-20g",$t);
+    printf(TEMP "%-20.2f",$t);
     for(my $k=0; ($k<$nrepl); $k++) {
        my $oo = shift;
        printf(TEMP "  %3d",$oo);

diff --git a/src/gromacs/gmxana/gmx_disre.cpp b/src/gromacs/gmxana/gmx_disre.cpp
index 9c6fe6564fbaf32d275c7f4d015389c5466a302e..2bed1d997f473640a614cbe56a00f5b2d839e9ff 100644 (file)
--- a/src/gromacs/gmxana/gmx_disre.cpp
+++ b/src/gromacs/gmxana/gmx_disre.cpp
@@ -41,6 +41,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <unordered_map>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/commandline/viewit.h"
@@ -181,10 +182,6 @@ static void check_viol(FILE *log,
             fprintf(debug, "DISRE: ndr = %d, label = %d  i=%d, n =%d\n",
                     ndr, label, i, n);
         }
-        if (ndr != label)
-        {
-            gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
-        }
         do
         {
             n += nat;
@@ -256,7 +253,7 @@ static void check_viol(FILE *log,
 }
 
 typedef struct {
-    int      index;
+    int      label;
     gmx_bool bCore;
     real     up1, r, rT3, rT6, viol, violT3, violT6;
 } t_dr_stats;
@@ -335,7 +332,7 @@ static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
             break;
         }
         fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
-                drs[i].index, yesno_names[drs[i].bCore],
+                drs[i].label, yesno_names[drs[i].bCore],
                 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
                 drs[i].viol, drs[i].violT3, drs[i].violT6);
     }
@@ -354,30 +351,26 @@ static gmx_bool is_core(int i, int isize, const int index[])
     return bIC;
 }
 
-static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
+static void dump_stats(FILE *log, int nsteps,
+                       const t_disresdata &dd,
+                       const t_ilist *disres,
                        t_iparams ip[], t_dr_result *dr,
                        int isize, int index[], t_atoms *atoms)
 {
-    int         i, j, nra;
     t_dr_stats *drs;
 
     fprintf(log, "\n");
     fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
-    snew(drs, ndr);
-    j         = 0;
-    nra       = interaction_function[F_DISRES].nratoms+1;
-    for (i = 0; (i < ndr); i++)
+    snew(drs, dd.nres);
+    const int nra = interaction_function[F_DISRES].nratoms + 1;
+    for (int j = 0; j < disres->nr; j += nra)
     {
-        /* Search for the appropriate restraint */
-        for (; (j < disres->nr); j += nra)
-        {
-            if (ip[disres->iatoms[j]].disres.label == i)
-            {
-                break;
-            }
-        }
-        drs[i].index  = i;
-        drs[i].bCore  = is_core(i, isize, index);
+        // Note that the restraint i can be used by multiple pairs
+        const int i   = disres->iatoms[j] - dd.type_min;
+        GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
+
+        drs[i].label  = ip[disres->iatoms[j]].disres.label;
+        drs[i].bCore  = is_core(drs[i].label, isize, index);
         drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
         drs[i].r      = dr->aver1[i]/nsteps;
         drs[i].rT3    = gmx::invcbrt(dr->aver_3[i]/nsteps);
@@ -393,23 +386,25 @@ static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
             atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
         }
     }
-    dump_viol(log, ndr, drs, FALSE);
+    dump_viol(log, dd.nres, drs, FALSE);
 
     fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
-    std::sort(drs, drs+ndr, [](const t_dr_stats &a, const t_dr_stats &b)
+    std::sort(drs, drs + dd.nres, [](const t_dr_stats &a, const t_dr_stats &b)
               {return a.viol > b.viol; });            //Reverse sort
-    dump_viol(log, ndr, drs, TRUE);
+    dump_viol(log, dd.nres, drs, TRUE);
 
-    dump_dump(log, ndr, drs);
+    dump_dump(log, dd.nres, drs);
 
     sfree(drs);
 }
 
-static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
+static void dump_clust_stats(FILE *fp,
+                             const t_disresdata &dd,
+                             const t_ilist *disres,
                              t_iparams ip[], t_blocka *clust, t_dr_result dr[],
                              char *clust_name[], int isize, int index[])
 {
-    int         i, j, k, nra, mmm = 0;
+    int         k, nra, mmm = 0;
     double      sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
     t_dr_stats *drs;
 
@@ -417,7 +412,7 @@ static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
     fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
     fprintf(fp, "Cluster  NFrames    SumV      MaxV     SumVT     MaxVT     SumVS     MaxVS\n");
 
-    snew(drs, ndr);
+    snew(drs, dd.nres);
 
     for (k = 0; (k < clust->nr); k++)
     {
@@ -436,21 +431,23 @@ static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
         {
             gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
         }
-        j         = 0;
         nra       = interaction_function[F_DISRES].nratoms+1;
         sumV      = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
-        for (i = 0; (i < ndr); i++)
+
+        // Use a map to process each restraint only once while looping over all pairs
+        std::unordered_map<int, bool> restraintHasBeenProcessed;
+        for (int j = 0; j < dd.nres; j += nra)
         {
-            /* Search for the appropriate restraint */
-            for (; (j < disres->nr); j += nra)
+            // Note that the restraint i can be used by multiple pairs
+            const int i   = disres->iatoms[j] - dd.type_min;
+
+            if (restraintHasBeenProcessed[i])
             {
-                if (ip[disres->iatoms[j]].disres.label == i)
-                {
-                    break;
-                }
+                continue;
             }
-            drs[i].index  = i;
-            drs[i].bCore  = is_core(i, isize, index);
+
+            drs[i].label  = ip[disres->iatoms[j]].disres.label;
+            drs[i].bCore  = is_core(drs[i].label, isize, index);
             drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
             drs[i].r      = dr[k].aver1[i]/dr[k].nframes;
             if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
@@ -468,6 +465,9 @@ static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
             maxV          = std::max(maxV, static_cast<double>(drs[i].viol));
             maxVT3        = std::max(maxVT3, static_cast<double>(drs[i].violT3));
             maxVT6        = std::max(maxVT6, static_cast<double>(drs[i].violT6));
+
+            // We have processed restraint i, mark it as such
+            restraintHasBeenProcessed[i] = true;
         }
         if (std::strcmp(clust_name[k], "1000") == 0)
         {
@@ -631,7 +631,7 @@ int gmx_disre(int argc, char *argv[])
         "because this analysis is done from a trajectory file afterwards",
         "it does not make sense to use time averaging. However,",
         "the time averaged values per restraint are given in the log file.[PAR]",
-        "An index file may be used to select specific restraints for",
+        "An index file may be used to select specific restraints by index group label for",
         "printing.[PAR]",
         "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
         "amount of average violations.[PAR]",
@@ -886,13 +886,13 @@ int gmx_disre(int argc, char *argv[])
 
     if (clust)
     {
-        dump_clust_stats(fplog, fcd.disres.nres, &(top.idef.il[F_DISRES]),
+        dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]),
                          top.idef.iparams, clust->clust, dr_clust,
                          clust->grpname, isize, index);
     }
     else
     {
-        dump_stats(fplog, j, fcd.disres.nres, &(top.idef.il[F_DISRES]),
+        dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]),
                    top.idef.iparams, &dr, isize, index,
                    bPDB ? atoms.get() : nullptr);
         if (bPDB)

diff --git a/src/gromacs/gpu_utils/cudautils.cuh b/src/gromacs/gpu_utils/cudautils.cuh
index 35ab3d875346b302c0926a938bee155d3976aab6..3efb65d0a9e8db55ef218242194173b1b9efdd2e 100644 (file)
--- a/src/gromacs/gpu_utils/cudautils.cuh
+++ b/src/gromacs/gpu_utils/cudautils.cuh
@@ -304,13 +304,14 @@ void prepareGpuKernelArgument(KernelPtr kernel,
  * A wrapper function for setting up all the CUDA kernel arguments.
  * Calls the recursive functions above.
  *
+ * \tparam    KernelPtr       Kernel function handle type
  * \tparam    Args            Types of all the kernel arguments
  * \param[in] kernel          Kernel function handle
  * \param[in] argsPtrs        Pointers to all the kernel arguments
  * \returns A prepared parameter pack to be used with launchGpuKernel() as the last argument.
  */
-template <typename ... Args>
-std::array<void *, sizeof ... (Args)> prepareGpuKernelArguments(void                     (*kernel)(Args...),
+template <typename KernelPtr, typename ... Args>
+std::array<void *, sizeof ... (Args)> prepareGpuKernelArguments(KernelPtr                kernel,
                                                                 const KernelLaunchConfig & /*config */,
                                                                 const Args *...          argsPtrs)
 {