Merge branch release-5-1 into release-2016
[alexxy/gromacs.git] / src / gromacs / domdec / domdec_topology.cpp
index 53b86035ec2dee927826cb5c23138d4f0db6bae0..46b05ddd71ddf01e40ece4a02723cd8fa42f0143 100644 (file)
 #include <string.h>
 
 #include <algorithm>
+#include <string>
 
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_network.h"
-#include "gromacs/legacyheaders/chargegroup.h"
-#include "gromacs/legacyheaders/force.h"
-#include "gromacs/legacyheaders/gmx_ga2la.h"
-#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/network.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/vsite.h"
-#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/domdec/ga2la.h"
+#include "gromacs/gmxlib/chargegroup.h"
+#include "gromacs/gmxlib/network.h"
 #include "gromacs/math/vec.h"
+#include "gromacs/mdlib/force.h"
+#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/mshift.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/topology/topsort.h"
 #include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
 
 #include "domdec_constraints.h"
 #include "domdec_internal.h"
@@ -102,7 +107,8 @@ typedef struct {
 } thread_work_t;
 
 /*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
-typedef struct gmx_reverse_top {
+struct gmx_reverse_top_t
+{
     //! @cond Doxygen_Suppress
     gmx_bool         bExclRequired;               /**< Do we require all exclusions to be assigned? */
     int              n_excl_at_max;               /**< The maximum number of exclusions one atom can have */
@@ -122,12 +128,8 @@ typedef struct gmx_reverse_top {
     /* Work data structures for multi-threading */
     int            nthread;           /**< The number of threads to be used */
     thread_work_t *th_work;           /**< Thread work array for local topology generation */
-
-    /* Pointers only used for an error message */
-    gmx_mtop_t     *err_top_global; /**< Pointer to the global top, only used for error reporting */
-    gmx_localtop_t *err_top_local;  /**< Pointer to the local top, only used for error reporting */
     //! @endcond
-} gmx_reverse_top_t;
+};
 
 /*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */
 static int nral_rt(int ftype)
@@ -158,7 +160,7 @@ static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck,
 }
 
 /*! \brief Print a header on error messages */
-static void print_error_header(FILE *fplog, char *moltypename, int nprint)
+static void print_error_header(FILE *fplog, const char *moltypename, int nprint)
 {
     fprintf(fplog, "\nMolecule type '%s'\n", moltypename);
     fprintf(stderr, "\nMolecule type '%s'\n", moltypename);
@@ -172,47 +174,41 @@ static void print_error_header(FILE *fplog, char *moltypename, int nprint)
 
 /*! \brief Help print error output when interactions are missing */
 static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
-                                          gmx_reverse_top_t *rt,
-                                          char *moltypename,
-                                          reverse_ilist_t *ril,
+                                          const gmx_reverse_top_t *rt,
+                                          const char *moltypename,
+                                          const reverse_ilist_t *ril,
                                           int a_start, int a_end,
                                           int nat_mol, int nmol,
-                                          t_idef *idef)
+                                          const t_idef *idef)
 {
-    int      nril_mol, *assigned, *gatindex;
-    int      ftype, ftype_j, nral, i, j_mol, j, a0, a0_mol, mol, a;
-    int      nprint;
-    t_ilist *il;
-    t_iatom *ia;
-    gmx_bool bFound;
-
-    nril_mol = ril->index[nat_mol];
+    int *assigned;
+    int  nril_mol = ril->index[nat_mol];
     snew(assigned, nmol*nril_mol);
 
-    gatindex = cr->dd->gatindex;
-    for (ftype = 0; ftype < F_NRE; ftype++)
+    int *gatindex = cr->dd->gatindex;
+    for (int ftype = 0; ftype < F_NRE; ftype++)
     {
         if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle))
         {
-            nral = NRAL(ftype);
-            il   = &idef->il[ftype];
-            ia   = il->iatoms;
-            for (i = 0; i < il->nr; i += 1+nral)
+            int            nral = NRAL(ftype);
+            const t_ilist *il   = &idef->il[ftype];
+            const t_iatom *ia   = il->iatoms;
+            for (int i = 0; i < il->nr; i += 1+nral)
             {
-                a0     = gatindex[ia[1]];
+                int a0 = gatindex[ia[1]];
                 /* Check if this interaction is in
                  * the currently checked molblock.
                  */
                 if (a0 >= a_start && a0 < a_end)
                 {
-                    mol    = (a0 - a_start)/nat_mol;
-                    a0_mol = (a0 - a_start) - mol*nat_mol;
-                    j_mol  = ril->index[a0_mol];
-                    bFound = FALSE;
-                    while (j_mol < ril->index[a0_mol+1] && !bFound)
+                    int  mol    = (a0 - a_start)/nat_mol;
+                    int  a0_mol = (a0 - a_start) - mol*nat_mol;
+                    int  j_mol  = ril->index[a0_mol];
+                    bool found  = false;
+                    while (j_mol < ril->index[a0_mol+1] && !found)
                     {
-                        j       = mol*nril_mol + j_mol;
-                        ftype_j = ril->il[j_mol];
+                        int j       = mol*nril_mol + j_mol;
+                        int ftype_j = ril->il[j_mol];
                         /* Here we need to check if this interaction has
                          * not already been assigned, since we could have
                          * multiply defined interactions.
@@ -221,23 +217,23 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
                             assigned[j] == 0)
                         {
                             /* Check the atoms */
-                            bFound = TRUE;
-                            for (a = 0; a < nral; a++)
+                            found = true;
+                            for (int a = 0; a < nral; a++)
                             {
                                 if (gatindex[ia[1+a]] !=
                                     a_start + mol*nat_mol + ril->il[j_mol+2+a])
                                 {
-                                    bFound = FALSE;
+                                    found = false;
                                 }
                             }
-                            if (bFound)
+                            if (found)
                             {
                                 assigned[j] = 1;
                             }
                         }
                         j_mol += 2 + nral_rt(ftype_j);
                     }
-                    if (!bFound)
+                    if (!found)
                     {
                         gmx_incons("Some interactions seem to be assigned multiple times");
                     }
@@ -249,16 +245,16 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
 
     gmx_sumi(nmol*nril_mol, assigned, cr);
 
-    nprint = 10;
-    i      = 0;
-    for (mol = 0; mol < nmol; mol++)
+    int nprint = 10;
+    int i      = 0;
+    for (int mol = 0; mol < nmol; mol++)
     {
-        j_mol = 0;
+        int j_mol = 0;
         while (j_mol < nril_mol)
         {
-            ftype = ril->il[j_mol];
-            nral  = NRAL(ftype);
-            j     = mol*nril_mol + j_mol;
+            int ftype = ril->il[j_mol];
+            int nral  = NRAL(ftype);
+            int j     = mol*nril_mol + j_mol;
             if (assigned[j] == 0 &&
                 !(interaction_function[ftype].flags & IF_VSITE))
             {
@@ -272,6 +268,7 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
                             interaction_function[ftype].longname);
                     fprintf(stderr, "%20s atoms",
                             interaction_function[ftype].longname);
+                    int a;
                     for (a = 0; a < nral; a++)
                     {
                         fprintf(fplog, "%5d", ril->il[j_mol+2+a]+1);
@@ -310,11 +307,12 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
 
 /*! \brief Help print error output when interactions are missing */
 static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
-                                             gmx_mtop_t *mtop, t_idef *idef)
+                                             const gmx_mtop_t *mtop,
+                                             const t_idef *idef)
 {
-    int                mb, a_start, a_end;
-    gmx_molblock_t    *molb;
-    gmx_reverse_top_t *rt;
+    int                      mb, a_start, a_end;
+    const gmx_molblock_t    *molb;
+    const gmx_reverse_top_t *rt;
 
     rt = cr->dd->reverse_top;
 
@@ -335,20 +333,19 @@ static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
     }
 }
 
-void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,  gmx_mtop_t *top_global, t_state *state_local)
+void dd_print_missing_interactions(FILE *fplog, t_commrec *cr,
+                                   int local_count,
+                                   const gmx_mtop_t *top_global,
+                                   const gmx_localtop_t *top_local,
+                                   t_state *state_local)
 {
     int             ndiff_tot, cl[F_NRE], n, ndiff, rest_global, rest_local;
     int             ftype, nral;
     char            buf[STRLEN];
     gmx_domdec_t   *dd;
-    gmx_mtop_t     *err_top_global;
-    gmx_localtop_t *err_top_local;
 
     dd = cr->dd;
 
-    err_top_global = dd->reverse_top->err_top_global;
-    err_top_local  = dd->reverse_top->err_top_local;
-
     if (fplog)
     {
         fprintf(fplog, "\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
@@ -360,7 +357,7 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
     for (ftype = 0; ftype < F_NRE; ftype++)
     {
         nral      = NRAL(ftype);
-        cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
+        cl[ftype] = top_local->idef.il[ftype].nr/(1+nral);
     }
 
     gmx_sumi(F_NRE, cl, cr);
@@ -386,10 +383,10 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
                 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
                 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
             {
-                n    = gmx_mtop_ftype_count(err_top_global, ftype);
+                n    = gmx_mtop_ftype_count(top_global, ftype);
                 if (ftype == F_CONSTR)
                 {
-                    n += gmx_mtop_ftype_count(err_top_global, F_CONSTRNC);
+                    n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
                 }
                 ndiff = cl[ftype] - n;
                 if (ndiff != 0)
@@ -420,21 +417,21 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
         }
     }
 
-    print_missing_interactions_atoms(fplog, cr, err_top_global,
-                                     &err_top_local->idef);
+    print_missing_interactions_atoms(fplog, cr, top_global, &top_local->idef);
     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
                  -1, state_local->x, state_local->box);
-    if (DDMASTER(dd))
+
+    std::string errorMessage;
+
+    if (ndiff_tot > 0)
     {
-        if (ndiff_tot > 0)
-        {
-            gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
-        }
-        else
-        {
-            gmx_fatal(FARGS, "%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_multibody(cr->dd), dd_cutoff_twobody(cr->dd));
-        }
+        errorMessage = "One or more interactions were assigned to multiple domains of the domain decompostion. Please report this bug.";
+    }
+    else
+    {
+        errorMessage = gmx::formatString("%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_multibody(dd), dd_cutoff_twobody(dd));
     }
+    gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), errorMessage.c_str());
 }
 
 /*! \brief Return global topology molecule information for global atom index \p i_gl */
@@ -507,9 +504,8 @@ static void count_excls(const t_block *cgs, const t_blocka *excls,
 }
 
 /*! \brief Run the reverse ilist generation and store it when \p bAssign = TRUE */
-static int low_make_reverse_ilist(const t_ilist *il_mt,
-                                  const t_atom *atom,
-                                  int **vsite_pbc, /* should be const */
+static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom,
+                                  const int * const * vsite_pbc,
                                   int *count,
                                   gmx_bool bConstr, gmx_bool bSettle,
                                   gmx_bool bBCheck,
@@ -519,8 +515,8 @@ static int low_make_reverse_ilist(const t_ilist *il_mt,
 {
     int            ftype, nral, i, j, nlink, link;
     const t_ilist *il;
-    t_iatom       *ia;
-    atom_id        a;
+    const t_iatom *ia;
+    int            a;
     int            nint;
     gmx_bool       bVSite;
 
@@ -613,7 +609,7 @@ static int low_make_reverse_ilist(const t_ilist *il_mt,
 /*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
 static int make_reverse_ilist(const t_ilist *ilist,
                               const t_atoms *atoms,
-                              int **vsite_pbc, /* should be const (C issue) */
+                              const int * const * vsite_pbc,
                               gmx_bool bConstr, gmx_bool bSettle,
                               gmx_bool bBCheck,
                               gmx_bool bLinkToAllAtoms,
@@ -659,8 +655,8 @@ static void destroy_reverse_ilist(reverse_ilist_t *ril)
 }
 
 /*! \brief Generate the reverse topology */
-static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE,
-                                           int ***vsite_pbc_molt,
+static gmx_reverse_top_t *make_reverse_top(const gmx_mtop_t *mtop, gmx_bool bFE,
+                                           const int * const * const * vsite_pbc_molt,
                                            gmx_bool bConstr, gmx_bool bSettle,
                                            gmx_bool bBCheck, int *nint)
 {
@@ -766,9 +762,9 @@ static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE,
 }
 
 void dd_make_reverse_top(FILE *fplog,
-                         gmx_domdec_t *dd, gmx_mtop_t *mtop,
-                         gmx_vsite_t *vsite,
-                         t_inputrec *ir, gmx_bool bBCheck)
+                         gmx_domdec_t *dd, const gmx_mtop_t *mtop,
+                         const gmx_vsite_t *vsite,
+                         const t_inputrec *ir, gmx_bool bBCheck)
 {
     if (fplog)
     {
@@ -813,7 +809,7 @@ void dd_make_reverse_top(FILE *fplog,
      * excluded pair should appear exactly once.
      */
     rt->bExclRequired = (ir->cutoff_scheme == ecutsGROUP &&
-                         IR_EXCL_FORCES(*ir));
+                         inputrecExclForces(ir));
 
     int nexcl, mb;
 
@@ -875,7 +871,7 @@ void dd_make_reverse_top(FILE *fplog,
  * atom-indexing organization code with the ifunc-adding code, so that
  * they can see that nral is the same value. */
 static gmx_inline void
-add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t ga2la,
+add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t *ga2la,
                      int nral, gmx_bool bHomeA,
                      int a, int a_gl, int a_mol,
                      const t_iatom *iatoms,
@@ -1026,7 +1022,7 @@ static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
 }
 
 /*! \brief Store a virtual site interaction, complex because of PBC and recursion */
-static void add_vsite(gmx_ga2la_t ga2la, const int *index, const int *rtil,
+static void add_vsite(gmx_ga2la_t *ga2la, const int *index, const int *rtil,
                       int ftype, int nral,
                       gmx_bool bHomeA, int a, int a_gl, int a_mol,
                       const t_iatom *iatoms,
@@ -1661,11 +1657,11 @@ static int make_exclusions_zone_cg(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                    int iz,
                                    int cg_start, int cg_end)
 {
-    int             n_excl_at_max, n, count, jla0, jla1, jla;
-    int             cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
-    const t_blocka *excls;
-    gmx_ga2la_t     ga2la;
-    int             cell;
+    int               n_excl_at_max, n, count, jla0, jla1, jla;
+    int               cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
+    const t_blocka   *excls;
+    gmx_ga2la_t      *ga2la;
+    int               cell;
 
     ga2la = dd->ga2la;
 
@@ -1798,9 +1794,9 @@ static void make_exclusions_zone(gmx_domdec_t *dd,
                                  int iz,
                                  int at_start, int at_end)
 {
-    gmx_ga2la_t ga2la;
-    int         jla0, jla1;
-    int         n_excl_at_max, n, at;
+    gmx_ga2la_t *ga2la;
+    int          jla0, jla1;
+    int          n_excl_at_max, n, at;
 
     ga2la = dd->ga2la;
 
@@ -1996,88 +1992,92 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd,
 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
         for (thread = 0; thread < rt->nthread; thread++)
         {
-            int       cg0t, cg1t;
-            t_idef   *idef_t;
-            int     **vsite_pbc;
-            int      *vsite_pbc_nalloc;
-            t_blocka *excl_t;
-
-            cg0t = cg0 + ((cg1 - cg0)* thread   )/rt->nthread;
-            cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
-
-            if (thread == 0)
+            try
             {
-                idef_t = idef;
-            }
-            else
-            {
-                idef_t = &rt->th_work[thread].idef;
-                clear_idef(idef_t);
-            }
+                int       cg0t, cg1t;
+                t_idef   *idef_t;
+                int     **vsite_pbc;
+                int      *vsite_pbc_nalloc;
+                t_blocka *excl_t;
+
+                cg0t = cg0 + ((cg1 - cg0)* thread   )/rt->nthread;
+                cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
 
-            if (vsite && vsite->bHaveChargeGroups && vsite->n_intercg_vsite > 0)
-            {
                 if (thread == 0)
                 {
-                    vsite_pbc        = vsite->vsite_pbc_loc;
-                    vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
+                    idef_t = idef;
                 }
                 else
                 {
-                    vsite_pbc        = rt->th_work[thread].vsite_pbc;
-                    vsite_pbc_nalloc = rt->th_work[thread].vsite_pbc_nalloc;
+                    idef_t = &rt->th_work[thread].idef;
+                    clear_idef(idef_t);
                 }
-            }
-            else
-            {
-                vsite_pbc        = NULL;
-                vsite_pbc_nalloc = NULL;
-            }
 
-            rt->th_work[thread].nbonded =
-                make_bondeds_zone(dd, zones,
-                                  mtop->molblock,
-                                  bRCheckMB, rcheck, bRCheck2B, rc2,
-                                  la2lc, pbc_null, cg_cm, idef->iparams,
-                                  idef_t,
-                                  vsite_pbc, vsite_pbc_nalloc,
-                                  izone,
-                                  dd->cgindex[cg0t], dd->cgindex[cg1t]);
-
-            if (izone < nzone_excl)
-            {
-                if (thread == 0)
+                if (vsite && vsite->bHaveChargeGroups && vsite->n_intercg_vsite > 0)
                 {
-                    excl_t = lexcls;
+                    if (thread == 0)
+                    {
+                        vsite_pbc        = vsite->vsite_pbc_loc;
+                        vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
+                    }
+                    else
+                    {
+                        vsite_pbc        = rt->th_work[thread].vsite_pbc;
+                        vsite_pbc_nalloc = rt->th_work[thread].vsite_pbc_nalloc;
+                    }
                 }
                 else
                 {
-                    excl_t      = &rt->th_work[thread].excl;
-                    excl_t->nr  = 0;
-                    excl_t->nra = 0;
+                    vsite_pbc        = NULL;
+                    vsite_pbc_nalloc = NULL;
                 }
 
-                if (dd->cgindex[dd->ncg_tot] == dd->ncg_tot &&
-                    !rt->bExclRequired)
-                {
-                    /* No charge groups and no distance check required */
-                    make_exclusions_zone(dd, zones,
-                                         mtop->moltype, cginfo,
-                                         excl_t,
-                                         izone,
-                                         cg0t, cg1t);
-                }
-                else
+                rt->th_work[thread].nbonded =
+                    make_bondeds_zone(dd, zones,
+                                      mtop->molblock,
+                                      bRCheckMB, rcheck, bRCheck2B, rc2,
+                                      la2lc, pbc_null, cg_cm, idef->iparams,
+                                      idef_t,
+                                      vsite_pbc, vsite_pbc_nalloc,
+                                      izone,
+                                      dd->cgindex[cg0t], dd->cgindex[cg1t]);
+
+                if (izone < nzone_excl)
                 {
-                    rt->th_work[thread].excl_count =
-                        make_exclusions_zone_cg(dd, zones,
-                                                mtop->moltype, bRCheck2B, rc2,
-                                                la2lc, pbc_null, cg_cm, cginfo,
-                                                excl_t,
-                                                izone,
-                                                cg0t, cg1t);
+                    if (thread == 0)
+                    {
+                        excl_t = lexcls;
+                    }
+                    else
+                    {
+                        excl_t      = &rt->th_work[thread].excl;
+                        excl_t->nr  = 0;
+                        excl_t->nra = 0;
+                    }
+
+                    if (dd->cgindex[dd->ncg_tot] == dd->ncg_tot &&
+                        !rt->bExclRequired)
+                    {
+                        /* No charge groups and no distance check required */
+                        make_exclusions_zone(dd, zones,
+                                             mtop->moltype, cginfo,
+                                             excl_t,
+                                             izone,
+                                             cg0t, cg1t);
+                    }
+                    else
+                    {
+                        rt->th_work[thread].excl_count =
+                            make_exclusions_zone_cg(dd, zones,
+                                                    mtop->moltype, bRCheck2B, rc2,
+                                                    la2lc, pbc_null, cg_cm, cginfo,
+                                                    excl_t,
+                                                    izone,
+                                                    cg0t, cg1t);
+                    }
                 }
             }
+            GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
         }
 
         if (rt->nthread > 1)
@@ -2134,7 +2134,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                        t_forcerec *fr,
                        rvec *cgcm_or_x,
                        gmx_vsite_t *vsite,
-                       gmx_mtop_t *mtop, gmx_localtop_t *ltop)
+                       const gmx_mtop_t *mtop, gmx_localtop_t *ltop)
 {
     gmx_bool bRCheckMB, bRCheck2B;
     real     rc = -1;
@@ -2194,8 +2194,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
             make_la2lc(dd);
             if (fr->bMolPBC)
             {
-                set_pbc_dd(&pbc, fr->ePBC, dd, TRUE, box);
-                pbc_null = &pbc;
+                pbc_null = set_pbc_dd(&pbc, fr->ePBC, dd->nc, TRUE, box);
             }
             else
             {
@@ -2223,13 +2222,9 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
     }
 
     ltop->atomtypes  = mtop->atomtypes;
-
-    /* For an error message only */
-    dd->reverse_top->err_top_global = mtop;
-    dd->reverse_top->err_top_local  = ltop;
 }
 
-void dd_sort_local_top(gmx_domdec_t *dd, t_mdatoms *mdatoms,
+void dd_sort_local_top(gmx_domdec_t *dd, const t_mdatoms *mdatoms,
                        gmx_localtop_t *ltop)
 {
     if (dd->reverse_top->ilsort == ilsortNO_FE)
@@ -2242,7 +2237,7 @@ void dd_sort_local_top(gmx_domdec_t *dd, t_mdatoms *mdatoms,
     }
 }
 
-gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
+gmx_localtop_t *dd_init_local_top(const gmx_mtop_t *top_global)
 {
     gmx_localtop_t *top;
     int             i;
@@ -2332,7 +2327,7 @@ static int *make_at2cg(t_block *cgs)
     return at2cg;
 }
 
-t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd,
+t_blocka *make_charge_group_links(const gmx_mtop_t *mtop, gmx_domdec_t *dd,
                                   cginfo_mb_t *cginfo_mb)
 {
     gmx_bool            bExclRequired;
@@ -2597,7 +2592,7 @@ static void bonded_cg_distance_mol(gmx_moltype_t *molt, int *at2cg,
 /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
 static void bonded_distance_intermol(const t_ilist *ilists_intermol,
                                      gmx_bool bBCheck,
-                                     rvec *x, int ePBC, matrix box,
+                                     const rvec *x, int ePBC, matrix box,
                                      bonded_distance_t *bd_2b,
                                      bonded_distance_t *bd_mb)
 {
@@ -2643,10 +2638,11 @@ static void bonded_distance_intermol(const t_ilist *ilists_intermol,
 }
 
 //! Compute charge group centers of mass for molecule \p molt
-static void get_cgcm_mol(gmx_moltype_t *molt, gmx_ffparams_t *ffparams,
+static void get_cgcm_mol(const gmx_moltype_t *molt,
+                         const gmx_ffparams_t *ffparams,
                          int ePBC, t_graph *graph, matrix box,
-                         gmx_vsite_t *vsite,
-                         rvec *x, rvec *xs, rvec *cg_cm)
+                         const gmx_vsite_t *vsite,
+                         const rvec *x, rvec *xs, rvec *cg_cm)
 {
     int n, i;
 
@@ -2705,8 +2701,9 @@ static int have_vsite_molt(gmx_moltype_t *molt)
 }
 
 void dd_bonded_cg_distance(FILE *fplog,
-                           gmx_mtop_t *mtop,
-                           t_inputrec *ir, rvec *x, matrix box,
+                           const gmx_mtop_t *mtop,
+                           const t_inputrec *ir,
+                           const rvec *x, matrix box,
                            gmx_bool bBCheck,
                            real *r_2b, real *r_mb)
 {
@@ -2720,7 +2717,7 @@ void dd_bonded_cg_distance(FILE *fplog,
     bonded_distance_t  bd_2b = { 0, -1, -1, -1 };
     bonded_distance_t  bd_mb = { 0, -1, -1, -1 };
 
-    bExclRequired = IR_EXCL_FORCES(*ir);
+    bExclRequired = inputrecExclForces(ir);
 
     vsite = init_vsite(mtop, NULL, TRUE);