Merge branch release-5-1
[alexxy/gromacs.git] / src / gromacs / mdlib / constr.cpp
index 7a3a215039e9b3c3b95e755b1db79b6db4dabcc7..5b850d0f5673be171eeca00c67a3827b5cdcb5ff 100644 (file)
@@ -36,7 +36,7 @@
  */
 #include "gmxpre.h"
 
-#include "gromacs/legacyheaders/constr.h"
+#include "constr.h"
 
 #include <assert.h>
 #include <stdlib.h>
 #include <algorithm>
 
 #include "gromacs/domdec/domdec.h"
+#include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/essentialdynamics/edsam.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/pdbio.h"
-#include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/mdrun.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/nrnb.h"
-#include "gromacs/legacyheaders/splitter.h"
-#include "gromacs/legacyheaders/txtdump.h"
-#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
+#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/splitter.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pulling/pull.h"
 #include "gromacs/topology/block.h"
 #include "gromacs/topology/invblock.h"
 #include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/pleasecite.h"
 #include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/txtdump.h"
 
 typedef struct gmx_constr {
     int                ncon_tot;       /* The total number of constraints    */
@@ -92,18 +95,18 @@ typedef struct gmx_constr {
     tensor            *vir_r_m_dr_th; /* Thread local working data          */
     int               *settle_error;  /* Thread local working data          */
 
-    gmx_mtop_t        *warn_mtop;     /* Only used for printing warnings    */
+    const gmx_mtop_t  *warn_mtop;     /* Only used for printing warnings    */
 } t_gmx_constr;
 
 typedef struct {
-    atom_id iatom[3];
-    atom_id blocknr;
+    int iatom[3];
+    int blocknr;
 } t_sortblock;
 
 static int pcomp(const void *p1, const void *p2)
 {
     int          db;
-    atom_id      min1, min2, max1, max2;
+    int          min1, min2, max1, max2;
     t_sortblock *a1 = (t_sortblock *)p1;
     t_sortblock *a2 = (t_sortblock *)p2;
 
@@ -170,7 +173,7 @@ void too_many_constraint_warnings(int eConstrAlg, int warncount)
 }
 
 static void write_constr_pdb(const char *fn, const char *title,
-                             gmx_mtop_t *mtop,
+                             const gmx_mtop_t *mtop,
                              int start, int homenr, t_commrec *cr,
                              rvec x[], matrix box)
 {
@@ -225,7 +228,7 @@ static void write_constr_pdb(const char *fn, const char *title,
     gmx_fio_fclose(out);
 }
 
-static void dump_confs(FILE *fplog, gmx_int64_t step, gmx_mtop_t *mtop,
+static void dump_confs(FILE *fplog, gmx_int64_t step, const gmx_mtop_t *mtop,
                        int start, int homenr, t_commrec *cr,
                        rvec x[], rvec xprime[], matrix box)
 {
@@ -360,7 +363,9 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
          * by the constraint coordinate communication routine,
          * so that here we can use normal pbc.
          */
-        pbc_null = set_pbc_dd(&pbc, ir->ePBC, cr->dd, FALSE, box);
+        pbc_null = set_pbc_dd(&pbc, ir->ePBC,
+                              DOMAINDECOMP(cr) ? cr->dd->nc : nullptr,
+                              FALSE, box);
     }
     else
     {
@@ -458,28 +463,32 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 #pragma omp parallel for num_threads(nth) schedule(static)
                 for (th = 0; th < nth; th++)
                 {
-                    int start_th, end_th;
-
-                    if (th > 0)
+                    try
                     {
-                        clear_mat(constr->vir_r_m_dr_th[th]);
+                        int start_th, end_th;
 
-                        constr->settle_error[th] = -1;
-                    }
+                        if (th > 0)
+                        {
+                            clear_mat(constr->vir_r_m_dr_th[th]);
 
-                    start_th = (nsettle* th   )/nth;
-                    end_th   = (nsettle*(th+1))/nth;
-                    if (start_th >= 0 && end_th - start_th > 0)
-                    {
-                        csettle(constr->settled,
-                                end_th-start_th,
-                                settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
-                                pbc_null,
-                                x[0], xprime[0],
-                                invdt, v ? v[0] : NULL, calcvir_atom_end,
-                                th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
-                                th == 0 ? &settle_error : &constr->settle_error[th]);
+                            constr->settle_error[th] = -1;
+                        }
+
+                        start_th = (nsettle* th   )/nth;
+                        end_th   = (nsettle*(th+1))/nth;
+                        if (start_th >= 0 && end_th - start_th > 0)
+                        {
+                            csettle(constr->settled,
+                                    end_th-start_th,
+                                    settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
+                                    pbc_null,
+                                    x[0], xprime[0],
+                                    invdt, v ? v[0] : NULL, calcvir_atom_end,
+                                    th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
+                                    th == 0 ? &settle_error : &constr->settle_error[th]);
+                        }
                     }
+                    GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
                 }
                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
                 if (v != NULL)
@@ -498,26 +507,30 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 #pragma omp parallel for num_threads(nth) schedule(static)
                 for (th = 0; th < nth; th++)
                 {
-                    int start_th, end_th;
-
-                    if (th > 0)
+                    try
                     {
-                        clear_mat(constr->vir_r_m_dr_th[th]);
-                    }
-
-                    start_th = (nsettle* th   )/nth;
-                    end_th   = (nsettle*(th+1))/nth;
-
-                    if (start_th >= 0 && end_th - start_th > 0)
-                    {
-                        settle_proj(constr->settled, econq,
-                                    end_th-start_th,
-                                    settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
-                                    pbc_null,
-                                    x,
-                                    xprime, min_proj, calcvir_atom_end,
-                                    th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th]);
+                        int start_th, end_th;
+
+                        if (th > 0)
+                        {
+                            clear_mat(constr->vir_r_m_dr_th[th]);
+                        }
+
+                        start_th = (nsettle* th   )/nth;
+                        end_th   = (nsettle*(th+1))/nth;
+
+                        if (start_th >= 0 && end_th - start_th > 0)
+                        {
+                            settle_proj(constr->settled, econq,
+                                        end_th-start_th,
+                                        settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
+                                        pbc_null,
+                                        x,
+                                        xprime, min_proj, calcvir_atom_end,
+                                        th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th]);
+                        }
                     }
+                    GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
                 }
                 /* This is an overestimate */
                 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
@@ -641,11 +654,11 @@ real *constr_rmsd_data(struct gmx_constr *constr)
     }
 }
 
-real constr_rmsd(struct gmx_constr *constr, gmx_bool bSD2)
+real constr_rmsd(struct gmx_constr *constr)
 {
     if (constr->lincsd)
     {
-        return lincs_rmsd(constr->lincsd, bSD2);
+        return lincs_rmsd(constr->lincsd);
     }
     else
     {
@@ -654,14 +667,14 @@ real constr_rmsd(struct gmx_constr *constr, gmx_bool bSD2)
 }
 
 static void make_shake_sblock_serial(struct gmx_constr *constr,
-                                     t_idef *idef, t_mdatoms *md)
+                                     t_idef *idef, const t_mdatoms *md)
 {
     int          i, j, m, ncons;
     int          bstart, bnr;
     t_blocka     sblocks;
     t_sortblock *sb;
     t_iatom     *iatom;
-    atom_id     *inv_sblock;
+    int         *inv_sblock;
 
     /* Since we are processing the local topology,
      * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
@@ -767,8 +780,8 @@ static void make_shake_sblock_serial(struct gmx_constr *constr,
 }
 
 static void make_shake_sblock_dd(struct gmx_constr *constr,
-                                 t_ilist *ilcon, t_block *cgs,
-                                 gmx_domdec_t *dd)
+                                 const t_ilist *ilcon, const t_block *cgs,
+                                 const gmx_domdec_t *dd)
 {
     int      ncons, c, cg;
     t_iatom *iatom;
@@ -799,7 +812,7 @@ static void make_shake_sblock_dd(struct gmx_constr *constr,
 }
 
 t_blocka make_at2con(int start, int natoms,
-                     t_ilist *ilist, t_iparams *iparams,
+                     const t_ilist *ilist, const t_iparams *iparams,
                      gmx_bool bDynamics, int *nflexiblecons)
 {
     int      *count, ncon, con, con_tot, nflexcon, ftype, i, a;
@@ -902,13 +915,13 @@ static int *make_at2settle(int natoms, const t_ilist *ilist)
 }
 
 void set_constraints(struct gmx_constr *constr,
-                     gmx_localtop_t *top, t_inputrec *ir,
-                     t_mdatoms *md, t_commrec *cr)
+                     gmx_localtop_t *top, const t_inputrec *ir,
+                     const t_mdatoms *md, t_commrec *cr)
 {
-    t_idef  *idef;
-    int      ncons;
-    t_ilist *settle;
-    int      iO, iH;
+    t_idef        *idef;
+    int            ncons;
+    const t_ilist *settle;
+    int            iO, iH;
 
     idef = &top->idef;
 
@@ -963,8 +976,9 @@ void set_constraints(struct gmx_constr *constr,
     }
 }
 
-static void constr_recur(t_blocka *at2con,
-                         t_ilist *ilist, t_iparams *iparams, gmx_bool bTopB,
+static void constr_recur(const t_blocka *at2con,
+                         const t_ilist *ilist, const t_iparams *iparams,
+                         gmx_bool bTopB,
                          int at, int depth, int nc, int *path,
                          real r0, real r1, real *r2max,
                          int *count)
@@ -1058,8 +1072,9 @@ static void constr_recur(t_blocka *at2con,
     }
 }
 
-static real constr_r_max_moltype(gmx_moltype_t *molt, t_iparams *iparams,
-                                 t_inputrec *ir)
+static real constr_r_max_moltype(const gmx_moltype_t *molt,
+                                 const t_iparams     *iparams,
+                                 const t_inputrec    *ir)
 {
     int      natoms, nflexcon, *path, at, count;
 
@@ -1126,7 +1141,7 @@ static real constr_r_max_moltype(gmx_moltype_t *molt, t_iparams *iparams,
     return rmax;
 }
 
-real constr_r_max(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir)
+real constr_r_max(FILE *fplog, const gmx_mtop_t *mtop, const t_inputrec *ir)
 {
     int  mt;
     real rmax;
@@ -1148,7 +1163,7 @@ real constr_r_max(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir)
 }
 
 gmx_constr_t init_constraints(FILE *fplog,
-                              gmx_mtop_t *mtop, t_inputrec *ir,
+                              const gmx_mtop_t *mtop, const t_inputrec *ir,
                               gmx_edsam_t ed, t_state *state,
                               t_commrec *cr)
 {