Merge "Merge release-5-0 into master"
authorRoland Schulz <roland@utk.edu>
Wed, 14 May 2014 16:25:27 +0000 (18:25 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Wed, 14 May 2014 16:25:28 +0000 (18:25 +0200)
22 files changed:
cmake/gmxManageSharedLibraries.cmake
src/gromacs/fileio/libxdrf.c
src/gromacs/fileio/tngio.cpp
src/gromacs/gmxana/geminate.c
src/gromacs/gmxana/gmx_anaeig.c
src/gromacs/gmxana/gmx_current.c
src/gromacs/gmxana/gmx_mindist.c
src/gromacs/gmxana/gmx_rms.c
src/gromacs/gmxana/gmx_rmsdist.c
src/gromacs/gmxana/gmx_xpm2ps.c
src/gromacs/gmxlib/gmx_detect_hardware.c
src/gromacs/gmxlib/nonbonded/nb_free_energy.c
src/gromacs/gmxpreprocess/readir.c
src/gromacs/gmxpreprocess/resall.c
src/gromacs/mdlib/coupling.c
src/gromacs/mdlib/domdec.c
src/gromacs/mdlib/domdec_setup.c
src/gromacs/mdlib/nbnxn_search.c
src/gromacs/mdlib/pme.c
src/gromacs/mdlib/qmmm.c
src/gromacs/options/basicoptions.cpp
src/gromacs/utility/cstringutil.c

index cc73de31b447636ff325072a75fc7a9daabf4116..fced521e952c7c254a6ae553b181c617a8385992 100644 (file)
@@ -39,7 +39,7 @@
 ########################################################################
 # Determine the defaults (this block has no effect if the variables have
 # already been set)
-if((APPLE OR CYGWIN OR ${CMAKE_SYSTEM_NAME} MATCHES "Linux|.*BSD") AND NOT GMX_BUILD_MDRUN_ONLY)
+if((APPLE OR CYGWIN OR ${CMAKE_SYSTEM_NAME} MATCHES "Linux|.*BSD|GNU") AND NOT GMX_BUILD_MDRUN_ONLY)
     # Maybe Solaris should be here? Patch this if you know!
     SET(SHARED_LIBS_DEFAULT ON)
 elseif(WIN32 OR ${CMAKE_SYSTEM_NAME} MATCHES "BlueGene")
index 758c348ed6b9817331f7242cd380a69bc823c02c..eb520c2e1f75feb1db867ea8fa86332403d1d986 100644 (file)
@@ -1345,7 +1345,6 @@ xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
     return res;
 }
 
-
 int
 xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
 {
@@ -1383,7 +1382,7 @@ xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
         {
             return -1;
         }
-        if (fr != frame && abs(low-high) > header_size)
+        if (fr != frame && llabs(low-high) > header_size)
         {
             if (fr < frame)
             {
@@ -1517,10 +1516,9 @@ int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeek
            the current time and the target time is bigger than dt and above all the distance between high
            and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
            if we reached the solution */
-        if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
-                                                                              - time) >= dt && dt_sign >= 0)
-             || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
-                                                        > header_size))
+        if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) ||
+             ((t - time) >= dt && dt_sign >= 0) || ((time - t) >= -dt && dt_sign < 0)) &&
+            (llabs(low - high) > header_size))
         {
             if (dt >= 0 && dt_sign != -1)
             {
@@ -1558,7 +1556,7 @@ int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeek
         }
         else
         {
-            if (abs(low - high) <= header_size)
+            if (llabs(low - high) <= header_size)
             {
                 break;
             }
index 89620eb7daa354106fe517f5a79688e12588714f..761b97501980a2ccc677810dcae2089f7dde306b 100644 (file)
@@ -197,9 +197,9 @@ static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
 
     /* FIXME: The TNG atoms should contain mass and atomB info (for free
      * energy calculations), i.e. in when it's available in TNG (2.0). */
-    for (int atomIt = 0; atomIt < atoms->nr; atomIt++)
+    for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
     {
-        const t_atom *at = &atoms->atom[atomIt];
+        const t_atom *at = &atoms->atom[atomIndex];
         /* FIXME: Currently the TNG API can only add atoms belonging to a
          * residue and chain. Wait for TNG 2.0*/
         if (atoms->nres > 0)
@@ -226,7 +226,7 @@ static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
             {
                 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
             }
-            tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIt]), *(atoms->atomtype[atomIt]), &tngAtom);
+            tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
         }
     }
     tng_molecule_cnt_set(tng, *tngMol, numMolecules);
@@ -245,18 +245,18 @@ void gmx_tng_add_mtop(tng_trajectory_t  tng,
         return;
     }
 
-    for (int molIt = 0; molIt < mtop->nmolblock; molIt++)
+    for (int molIndex = 0; molIndex < mtop->nmolblock; molIndex++)
     {
         tng_molecule_t       tngMol  = NULL;
         const gmx_moltype_t *molType =
-            &mtop->moltype[mtop->molblock[molIt].type];
+            &mtop->moltype[mtop->molblock[molIndex].type];
 
         /* Add a molecule to the TNG trajectory with the same name as the
          * current molecule. */
         addTngMoleculeFromTopology(tng,
                                    *(molType->name),
                                    &molType->atoms,
-                                   mtop->molblock[molIt].nmol,
+                                   mtop->molblock[molIndex].nmol,
                                    &tngMol);
 
         /* Bonds have to be deduced from interactions (constraints etc). Different
@@ -488,8 +488,44 @@ void gmx_tng_prepare_md_writing(tng_trajectory_t  tng,
 }
 
 #ifdef GMX_USE_TNG
-/* Create a TNG molecule representing the selection groups
- * to write */
+/* Check if all atoms in the molecule system are selected
+ * by a selection group of type specified by gtype. */
+static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
+                                   const int         gtype)
+{
+    const gmx_moltype_t     *molType;
+    const t_atoms           *atoms;
+
+    /* Iterate through all atoms in the molecule system and
+     * check if they belong to a selection group of the
+     * requested type. */
+    for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
+    {
+        molType = &mtop->moltype[mtop->molblock[molIndex].type];
+
+        atoms = &molType->atoms;
+
+        for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
+        {
+            for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
+            {
+                if (ggrpnr(&mtop->groups, gtype, i) != 0)
+                {
+                    return FALSE;
+                }
+            }
+        }
+    }
+    return TRUE;
+}
+
+/* Create TNG molecules which will represent each of the selection
+ * groups for writing. But do that only if there is actually a
+ * specified selection group and it is not the whole system.
+ * TODO: Currently the only selection that is taken into account
+ * is egcCompressedX, but other selections should be added when
+ * e.g. writing energies is implemented.
+ */
 static void add_selection_groups(tng_trajectory_t  tng,
                                  const gmx_mtop_t *mtop)
 {
@@ -498,7 +534,7 @@ static void add_selection_groups(tng_trajectory_t  tng,
     const t_atom            *at;
     const t_resinfo         *resInfo;
     const t_ilist           *ilist;
-    int                      nAtoms      = 0, i = 0, j, molIt, atomIt, nameIndex;
+    int                      nameIndex;
     int                      atom_offset = 0;
     tng_molecule_t           mol, iterMol;
     tng_chain_t              chain;
@@ -508,6 +544,25 @@ static void add_selection_groups(tng_trajectory_t  tng,
     gmx_int64_t              nMols;
     char                    *groupName;
 
+    /* TODO: When the TNG molecules block is more flexible TNG selection
+     * groups should not need all atoms specified. It should be possible
+     * just to specify what molecules are selected (and/or which atoms in
+     * the molecule) and how many (if applicable). */
+
+    /* If no atoms are selected we do not need to create a
+     * TNG selection group molecule. */
+    if (mtop->groups.ngrpnr[egcCompressedX] == 0)
+    {
+        return;
+    }
+
+    /* If all atoms are selected we do not have to create a selection
+     * group molecule in the TNG molecule system. */
+    if (all_atoms_selected(mtop, egcCompressedX))
+    {
+        return;
+    }
+
     /* The name of the TNG molecule containing the selection group is the
      * same as the name of the selection group. */
     nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
@@ -516,16 +571,16 @@ static void add_selection_groups(tng_trajectory_t  tng,
     tng_molecule_alloc(tng, &mol);
     tng_molecule_name_set(tng, mol, groupName);
     tng_molecule_chain_add(tng, mol, "", &chain);
-    for (molIt = 0; molIt < mtop->nmoltype; molIt++)
+    for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
     {
-        molType = &mtop->moltype[mtop->molblock[molIt].type];
+        molType = &mtop->moltype[mtop->molblock[molIndex].type];
 
         atoms = &molType->atoms;
 
-        for (j = 0; j < mtop->molblock[molIt].nmol; j++)
+        for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
         {
             bool bAtomsAdded = FALSE;
-            for (atomIt = 0; atomIt < atoms->nr; atomIt++, i++)
+            for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
             {
                 char *res_name;
                 int   res_id;
@@ -534,7 +589,7 @@ static void add_selection_groups(tng_trajectory_t  tng,
                 {
                     continue;
                 }
-                at = &atoms->atom[atomIt];
+                at = &atoms->atom[atomIndex];
                 if (atoms->nres > 0)
                 {
                     resInfo = &atoms->resinfo[at->resind];
@@ -555,10 +610,9 @@ static void add_selection_groups(tng_trajectory_t  tng,
                      * original residue IDs - otherwise there might be conflicts. */
                     tng_chain_residue_add(tng, chain, res_name, &res);
                 }
-                tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIt]),
-                                          *(atoms->atomtype[atomIt]),
-                                          atom_offset + atomIt, &atom);
-                nAtoms++;
+                tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
+                                          *(atoms->atomtype[atomIndex]),
+                                          atom_offset + atomIndex, &atom);
                 bAtomsAdded = TRUE;
             }
             /* Add bonds. */
@@ -619,24 +673,17 @@ static void add_selection_groups(tng_trajectory_t  tng,
             atom_offset += atoms->nr;
         }
     }
-    if (nAtoms != i)
+    tng_molecule_existing_add(tng, &mol);
+    tng_molecule_cnt_set(tng, mol, 1);
+    tng_num_molecule_types_get(tng, &nMols);
+    for (gmx_int64_t k = 0; k < nMols; k++)
     {
-        tng_molecule_existing_add(tng, &mol);
-        tng_molecule_cnt_set(tng, mol, 1);
-        tng_num_molecule_types_get(tng, &nMols);
-        for (gmx_int64_t k = 0; k < nMols; k++)
+        tng_molecule_of_index_get(tng, k, &iterMol);
+        if (iterMol == mol)
         {
-            tng_molecule_of_index_get(tng, k, &iterMol);
-            if (iterMol == mol)
-            {
-                continue;
-            }
-            tng_molecule_cnt_set(tng, iterMol, 0);
+            continue;
         }
-    }
-    else
-    {
-        tng_molecule_free(tng, &mol);
+        tng_molecule_cnt_set(tng, iterMol, 0);
     }
 }
 #endif
@@ -645,7 +692,7 @@ void gmx_tng_set_compression_precision(tng_trajectory_t tng,
                                        real             prec)
 {
 #ifdef GMX_USE_TNG
-    tng_compression_precision_set(tng, 1.0/prec);
+    tng_compression_precision_set(tng, prec);
 #else
     GMX_UNUSED_VALUE(tng);
     GMX_UNUSED_VALUE(prec);
index 32d8f8c5a7861cd206a3b780f7973b3c8b5726af..4a0a998a8c419dfc3c73f04d84da888735d39d48 100644 (file)
@@ -718,10 +718,10 @@ extern void fixGemACF(double *ct, int len)
     bBad = FALSE;
 
     /* An acf of binary data must be one at t=0. */
-    if (abs(ct[0]-1.0) > 1e-6)
+    if (fabs(ct[0]-1.0) > 1e-6)
     {
         ct[0] = 1.0;
-        fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0));
+        fprintf(stderr, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct[0]-1.0));
     }
 
     for (i = 0; i < len; i++)
index 958b2719360fce76c89bf40017a8c0229c0d6c7b..014a58c147719b68717c9ffe4919883aa4a23235 100644 (file)
@@ -260,7 +260,7 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
         {
             for (i = 0; i < n; i++)
             {
-                if (bSplit && i > 0 && abs(x[i]) < 1e-5)
+                if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
                 {
                     if (output_env_get_print_xvgr_codes(oenv))
                     {
@@ -673,7 +673,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb
                            oenv);
         for (i = 0; i < nframes; i++)
         {
-            if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5)
+            if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
             {
                 fprintf(xvgrout, "&\n");
             }
@@ -761,7 +761,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb
             j = 0;
             for (i = 0; i < atoms.nr; i++)
             {
-                if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5)
+                if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
                 {
                     fprintf(out, "TER\n");
                     j = 0;
index 171a7a91c24a01075cead154acbab4a26a0306b1..5154ae2ece25295727c0c6eed6135497379e8e3c 100644 (file)
@@ -141,7 +141,7 @@ static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
         }
     }
 
-    if (abs(qall) > 0.01)
+    if (fabs(qall) > 0.01)
     {
         printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
         bNEU = FALSE;
index 33e54d3e134de781e6852dd03293d59c3369af3e..73833170002de10535ae28de81ab03b52b255c4b 100644 (file)
@@ -170,7 +170,7 @@ static void periodic_mindist_plot(const char *trxfn, const char *outfn,
             ind_mini = ind_min[0];
             ind_minj = ind_min[1];
         }
-        if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
+        if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
         {
             fprintf(out, "&\n");
         }
@@ -421,7 +421,7 @@ void dist_plot(const char *fn, const char *afile, const char *dfile,
     bFirst = TRUE;
     do
     {
-        if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
+        if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
         {
             fprintf(dist, "&\n");
             if (num)
index 28cd76faf3ae735397435e26cf1f19d5087ba850..46acd9435968232c0256f99e88275207efb8c886 100644 (file)
@@ -1137,7 +1137,7 @@ int gmx_rms(int argc, char *argv[])
     for (i = 0; (i < teller); i++)
     {
         if (bSplit && i > 0 &&
-            abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
+            fabs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
         {
             fprintf(fp, "&\n");
         }
@@ -1179,7 +1179,7 @@ int gmx_rms(int argc, char *argv[])
         }
         for (i = 0; (i < teller); i++)
         {
-            if (bSplit && i > 0 && abs(time[i]) < 1e-5)
+            if (bSplit && i > 0 && fabs(time[i]) < 1e-5)
             {
                 fprintf(fp, "&\n");
             }
index 2d6886d2399cb50b6d8a3c46cc24ebed4567e568..9e40d9c4c14ec01aab1996ec0eb2dee48c04e1da 100644 (file)
@@ -325,7 +325,7 @@ static int analyze_noe_equivalent(const char *eq_fn,
                     if (bEquiv)
                     {
                         /* set index for matching atom */
-                        noe_index[j] = groupnr;
+                        noe_index[i] = groupnr;
                         /* skip matching atom */
                         i = j;
                     }
@@ -343,7 +343,7 @@ static int analyze_noe_equivalent(const char *eq_fn,
                    This is supposed to cover all CH3 groups and the like */
                 anmi   = *atoms->atomname[index[i]];
                 anmil  = strlen(anmi);
-                bMatch = i < isize-3 && anmi[anmil-1] == '1';
+                bMatch = i <= isize-3 && anmi[anmil-1] == '1';
                 if (bMatch)
                 {
                     for (j = 1; j < 3; j++)
@@ -646,7 +646,7 @@ int gmx_rmsdist(int argc, char *argv[])
         "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
         "a set of equivalent atoms specified as residue number and name and",
         "atom name; e.g.:[PAR]",
-        "[TT]3 SER  HB1 3 SER  HB2[tt][PAR]",
+        "[TT]HB* 3 SER  HB1 3 SER  HB2[tt][PAR]",
         "Residue and atom names must exactly match those in the structure",
         "file, including case. Specifying non-sequential atoms is undefined."
 
@@ -787,7 +787,8 @@ int gmx_rmsdist(int argc, char *argv[])
     }
 
     /*do a first step*/
-    natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+    natom  = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+    teller = 0;
 
     do
     {
@@ -795,14 +796,13 @@ int gmx_rmsdist(int argc, char *argv[])
 
         rmsnow = rms_diff(isize, d, d_r);
         fprintf(fp, "%g  %g\n", t, rmsnow);
+        teller++;
     }
     while (read_next_x(oenv, status, &t, x, box));
     fprintf(stderr, "\n");
 
     gmx_ffclose(fp);
 
-    teller = nframes_read(status);
-
     close_trj(status);
 
     calc_rms(isize, teller, dtot, dtot2, mean, &meanmax, rms, &rmsmax, rmsc, &rmscmax);
index d230bca60f6bcb8c349ef60d0246c9f6c282d5e3..bc491707f86e6a5fe80ddbf2b9f1ba60c83f6075 100644 (file)
@@ -523,8 +523,8 @@ static void draw_zerolines(t_psdata out, real x0, real y0, real w,
                 xx = xx00+(x+0.7)*psr->xboxsize;
                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
                 if (x != 0 && x < mat[i].nx-1 &&
-                    abs(mat[i].axis_x[x]) <
-                    0.1*abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
+                    fabs(mat[i].axis_x[x]) <
+                    0.1*fabs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
                 {
                     ps_line (out, xx, yy00, xx, yy00+dy+2);
                 }
@@ -538,8 +538,8 @@ static void draw_zerolines(t_psdata out, real x0, real y0, real w,
                 yy = yy00+(y+0.7)*psr->yboxsize;
                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
                 if (y != 0 && y < mat[i].ny-1 &&
-                    abs(mat[i].axis_y[y]) <
-                    0.1*abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
+                    fabs(mat[i].axis_y[y]) <
+                    0.1*fabs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
                 {
                     ps_line (out, xx00, yy, xx00+w+2, yy);
                 }
@@ -1135,7 +1135,7 @@ void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
             }
             for (x = 0; x < mats[i].nx-1; x++)
             {
-                if (abs(mats[i].axis_x[x+1]) < 1e-5)
+                if (fabs(mats[i].axis_x[x+1]) < 1e-5)
                 {
                     for (y = 0; y < mats[i].ny; y++)
                     {
@@ -1145,7 +1145,7 @@ void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
             }
             for (y = 0; y < mats[i].ny-1; y++)
             {
-                if (abs(mats[i].axis_y[y+1]) < 1e-5)
+                if (fabs(mats[i].axis_y[y+1]) < 1e-5)
                 {
                     for (x = 0; x < mats[i].nx; x++)
                     {
index f3400d33818bb3ff6144addc0c2e9b3f77df7e0e..96ce3a1a9330ea2a4d490b7df13082469f38ede5 100644 (file)
@@ -582,7 +582,7 @@ static void gmx_detect_gpus(FILE *fplog, const t_commrec *cr)
 
         if (detect_cuda_gpus(&hwinfo_g->gpu_info, detection_error) != 0)
         {
-            if (detection_error != NULL && detection_error[0] != '\0')
+            if (detection_error[0] != '\0')
             {
                 sprintf(sbuf, ":\n      %s\n", detection_error);
             }
index 74cd5640efc725b56dec4d93b2deec0c1c2f7eb9..5414c668c6936a22a2238e2a079fba88c91e34ae 100644 (file)
@@ -795,6 +795,7 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
      * 12  flops per outer iteration
      * 150 flops per inner iteration
      */
+#pragma omp atomic
     inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150);
 }
 
index 8c3da5f6f07a3ffd7a8029c2dff0ad4891b361b6..8b1b8c823829f72d1a1facb7a344f7d85ee9ddcb 100644 (file)
@@ -1142,13 +1142,22 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         warning_note(wi, warn_buf);
     }
 
-    if (ir->coulombtype == eelPMESWITCH)
+    if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
     {
         if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
         {
-            sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
-                    eel_names[ir->coulombtype],
-                    ir->ewald_rtol);
+            real percentage  = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
+            sprintf(warn_buf, "The switching range for should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
+                    percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
+            warning(wi, warn_buf);
+        }
+    }
+
+    if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
+    {
+        if (ir->rvdw_switch == 0)
+        {
+            sprintf(warn_buf, "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones potential.  This suggests it was not set in the mdp, which can lead to large energy errors.  In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw switching range.");
             warning(wi, warn_buf);
         }
     }
index 01e77a8a06fc7d662df1744ac9a09082ee8c7349..47616b2483e46e713346bcf8fcd81a71d9495167 100644 (file)
@@ -79,20 +79,23 @@ gpp_atomtype_t read_atype(const char *ffdir, t_symtab *tab)
             /* Skip blank or comment-only lines */
             do
             {
-                fgets2(buf, STRLEN, in);
-                if (NULL != buf)
+                if (fgets2(buf, STRLEN, in) != NULL)
                 {
                     strip_comment(buf);
                     trim(buf);
                 }
             }
-            while (!feof(in) && NULL != buf && strlen(buf) == 0);
+            while (!feof(in) && strlen(buf) == 0);
 
-            if ((buf != NULL) && (sscanf(buf, "%s%lf", name, &m) == 2))
+            if (sscanf(buf, "%s%lf", name, &m) == 2)
             {
                 a->m = m;
                 add_atomtype(at, tab, a, name, nb, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
-                fprintf(stderr, "\rAtomtype %d", nratt+1);
+                fprintf(stderr, "\rAtomtype %d", ++nratt);
+            }
+            else
+            {
+                fprintf(stderr, "\nInvalid format: %s\n", buf);
             }
         }
         gmx_ffclose(in);
index 395f1d088e09f54f92be3e49ef904459007a6fa1..bfc6d636b9f5a18fb1f5191546d0a4820eca178d 100644 (file)
@@ -1355,58 +1355,43 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
     return ener_npt;
 }
 
-static real vrescale_gamdev(int ia,
+static real vrescale_gamdev(real ia,
                             gmx_int64_t step, gmx_int64_t *count,
                             gmx_int64_t seed1, gmx_int64_t seed2)
 /* Gamma distribution, adapted from numerical recipes */
 {
-    int    j;
     real   am, e, s, v1, v2, x, y;
     double rnd[2];
 
-    if (ia < 6)
+    assert(ia > 1);
+
+    do
     {
         do
         {
-            x = 1.0;
-            for (j = 1; j <= ia; j++)
+            do
             {
                 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
-                x *= rnd[0];
+                v1 = rnd[0];
+                v2 = 2.0*rnd[1] - 1.0;
             }
+            while (v1*v1 + v2*v2 > 1.0 ||
+                   v1*v1*GMX_REAL_MAX < 3.0*ia);
+            /* The last check above ensures that both x (3.0 > 2.0 in s)
+             * and the pre-factor for e do not go out of range.
+             */
+            y  = v2/v1;
+            am = ia - 1;
+            s  = sqrt(2.0*am + 1.0);
+            x  = s*y + am;
         }
-        while (x == 0);
-        x = -log(x);
-    }
-    else
-    {
-        do
-        {
-            do
-            {
-                do
-                {
-                    gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
-                    v1 = rnd[0];
-                    v2 = 2.0*rnd[1] - 1.0;
-                }
-                while (v1*v1 + v2*v2 > 1.0 ||
-                       v1*v1*GMX_REAL_MAX < 3.0*ia);
-                /* The last check above ensures that both x (3.0 > 2.0 in s)
-                 * and the pre-factor for e do not go out of range.
-                 */
-                y  = v2/v1;
-                am = ia - 1;
-                s  = sqrt(2.0*am + 1.0);
-                x  = s*y + am;
-            }
-            while (x <= 0.0);
-            e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
+        while (x <= 0.0);
 
-            gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
-        }
-        while (rnd[0] > e);
+        e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
+
+        gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
     }
+    while (rnd[0] > e);
 
     return x;
 }
@@ -1430,30 +1415,48 @@ static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
     return x*r;
 }
 
-static real vrescale_sumnoises(int nn,
+static real vrescale_sumnoises(real nn,
                                gmx_int64_t step, gmx_int64_t *count,
                                gmx_int64_t seed1, gmx_int64_t seed2)
 {
 /*
- * Returns the sum of n independent gaussian noises squared
+ * Returns the sum of nn independent gaussian noises squared
  * (i.e. equivalent to summing the square of the return values
- * of nn calls to gmx_rng_gaussian_real).xs
+ * of nn calls to gmx_rng_gaussian_real).
  */
-    real r, gauss;
-
-    r = 2.0*vrescale_gamdev(nn/2, step, count, seed1, seed2);
+    const real ndeg_tol = 0.0001;
+    real       r;
 
-    if (nn % 2 == 1)
+    if (nn < 2 + ndeg_tol)
     {
-        gauss = gaussian_count(step, count, seed1, seed2);
+        int  nn_int, i;
+        real gauss;
+
+        nn_int = (int)(nn + 0.5);
+
+        if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
+        {
+            gmx_fatal(FARGS, "The v-rescale thermostat was called with a group with #DOF=%f, but for #DOF<3 only integer #DOF are supported", nn + 1);
+        }
+
+        r = 0;
+        for (i = 0; i < nn_int; i++)
+        {
+            gauss = gaussian_count(step, count, seed1, seed2);
 
-        r += gauss*gauss;
+            r += gauss*gauss;
+        }
+    }
+    else
+    {
+        /* Use a gamma distribution for any real nn > 2 */
+        r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
     }
 
     return r;
 }
 
-static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut,
+static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
                                  gmx_int64_t step, gmx_int64_t seed)
 {
 /*
index e45c02926ca42ddbbd968e2814b4b5f482ae50c3..b44c0cb3cb7f6d1d4a6dafb4e0f016d22766e56c 100644 (file)
@@ -453,8 +453,14 @@ static const ivec dd_zp1[dd_zp1n] = {{0, 0, 2}};
 /* Factor to account for pressure scaling during nstlist steps */
 #define DD_PRES_SCALE_MARGIN 1.02
 
-/* Allowed performance loss before we DLB or warn */
-#define DD_PERF_LOSS 0.05
+/* Turn on DLB when the load imbalance causes this amount of total loss.
+ * There is a bit of overhead with DLB and it's difficult to achieve
+ * a load imbalance of less than 2% with DLB.
+ */
+#define DD_PERF_LOSS_DLB_ON  0.02
+
+/* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
+#define DD_PERF_LOSS_WARN    0.05
 
 #define DD_CELL_F_SIZE(dd, di) ((dd)->nc[(dd)->dim[(di)]]+1+(di)*2+1+(di))
 
@@ -5479,7 +5485,7 @@ static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
         fprintf(fplog, "\n");
         fprintf(stderr, "\n");
 
-        if (lossf >= DD_PERF_LOSS)
+        if (lossf >= DD_PERF_LOSS_WARN)
         {
             sprintf(buf,
                     "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
@@ -5495,7 +5501,7 @@ static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
             fprintf(fplog, "%s\n", buf);
             fprintf(stderr, "%s\n", buf);
         }
-        if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS)
+        if (npme > 0 && fabs(lossp) >= DD_PERF_LOSS_WARN)
         {
             sprintf(buf,
                     "NOTE: %.1f %% performance was lost because the PME nodes\n"
@@ -9381,7 +9387,7 @@ void dd_partition_system(FILE                *fplog,
                 if (DDMASTER(dd))
                 {
                     bTurnOnDLB =
-                        (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS);
+                        (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
                     if (debug)
                     {
                         fprintf(debug, "step %s, imb loss %f\n",
index a7bdaaca632f46dd4251047c8f1b99fcfbf444f0..df6bd8770f1fee3d24c1431cb8be03a6974eabf0 100644 (file)
@@ -719,7 +719,7 @@ real dd_choose_grid(FILE *fplog,
             cr->npmenodes = 0;
         }
 
-        if (cr->nnodes > 12)
+        if (nnodes_div > 12)
         {
             ldiv = largest_divisor(nnodes_div);
             /* Check if the largest divisor is more than nnodes^2/3 */
index 3570d120d4f7c60929baaf2eb1aa094b9a0b7235..946f9701bebad6270ee30001ac48f48bf638e3b3 100644 (file)
@@ -3308,8 +3308,11 @@ static void make_fep_list(const nbnxn_search_t    nbs,
     cj_ind_start = nbl_ci->cj_ind_start;
     cj_ind_end   = nbl_ci->cj_ind_end;
 
-    /* In worst case we have alternating energy groups and create npair lists */
-    nri_max = nbl->na_ci*(cj_ind_end - cj_ind_start);
+    /* In worst case we have alternating energy groups
+     * and create #atom-pair lists, which means we need the size
+     * of a cluster pair (na_ci*na_cj) times the number of cj's.
+     */
+    nri_max = nbl->na_ci*nbl->na_cj*(cj_ind_end - cj_ind_start);
     if (nlist->nri + nri_max > nlist->maxnri)
     {
         nlist->maxnri = over_alloc_large(nlist->nri + nri_max);
@@ -3445,6 +3448,7 @@ static void make_fep_list(const nbnxn_search_t    nbs,
 
             if (nlist->nrj > nlist->jindex[nri])
             {
+                /* Actually add this new, non-empty, list */
                 nlist->nri++;
                 nlist->jindex[nlist->nri] = nlist->nrj;
             }
@@ -3510,8 +3514,14 @@ static void make_fep_list_supersub(const nbnxn_search_t    nbs,
     cj4_ind_start = nbl_sci->cj4_ind_start;
     cj4_ind_end   = nbl_sci->cj4_ind_end;
 
-    /* No energy groups (yet), so we split lists in max_nrj_fep pairs */
-    nri_max = nbl->na_sc*(1 + ((cj4_ind_end - cj4_ind_start)*NBNXN_GPU_JGROUP_SIZE)/max_nrj_fep);
+    /* Here we process one super-cell, max #atoms na_sc, versus a list
+     * cj4 entries, each with max NBNXN_GPU_JGROUP_SIZE cj's, each
+     * of size na_cj atoms.
+     * On the GPU we don't support energy groups (yet).
+     * So for each of the na_sc i-atoms, we need max one FEP list
+     * for each max_nrj_fep j-atoms.
+     */
+    nri_max = nbl->na_sc*nbl->na_cj*(1 + ((cj4_ind_end - cj4_ind_start)*NBNXN_GPU_JGROUP_SIZE)/max_nrj_fep);
     if (nlist->nri + nri_max > nlist->maxnri)
     {
         nlist->maxnri = over_alloc_large(nlist->nri + nri_max);
@@ -3634,6 +3644,7 @@ static void make_fep_list_supersub(const nbnxn_search_t    nbs,
 
                 if (nlist->nrj > nlist->jindex[nri])
                 {
+                    /* Actually add this new, non-empty, list */
                     nlist->nri++;
                     nlist->jindex[nlist->nri] = nlist->nrj;
                 }
index d56b2a41e5254587f2d7025735666e1c73b74b50..6c7f1b981a8ed72620aac0a27b0208520fb19864 100644 (file)
@@ -5154,10 +5154,10 @@ int gmx_pme_do(gmx_pme_t pme,
             }
             if (flags & GMX_PME_SOLVE)
             {
-                int loop_count;
                 /* solve in k-space for our local cells */
 #pragma omp parallel num_threads(pme->nthread) private(thread)
                 {
+                    int loop_count;
                     thread = gmx_omp_get_thread_num();
                     if (thread == 0)
                     {
index a5582fab5839b1555860888d33870f274cc7eeff..4e85721b9a91e9c56616b52b91799254d599307e 100644 (file)
@@ -196,7 +196,7 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec
 #elif defined GMX_QMMM_GAUSSIAN
             QMener = call_gaussian(cr, fr, qm, mm, f, fshift);
 #elif defined GMX_QMMM_ORCA
-            QMener = call_orca(cr, fr, qm, mm, f, fshift);
+            QMener = call_orca(fr, qm, mm, f, fshift);
 #else
             gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 #endif
@@ -769,7 +769,7 @@ void init_QMMMrec(t_commrec  *cr,
 #elif defined GMX_QMMM_GAUSSIAN
             init_gaussian(cr, qr->qm[0], qr->mm);
 #elif defined GMX_QMMM_ORCA
-            init_orca(cr, qr->qm[0], qr->mm);
+            init_orca(qr->qm[0]);
 #else
             gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 #endif
index 1d1a5f1a74d034b8c3e27d050f4103870bca5a7f..f1af57aa32bed7ea3c075503fbfb457b610b6b42 100644 (file)
@@ -376,7 +376,7 @@ void FloatOptionStorage::convertValue(const std::string &value)
     double      dval = std::strtod(ptr, &endptr);
     if (errno == ERANGE
         || dval * factor_ < -std::numeric_limits<float>::max()
-        || dval * factor_ > -std::numeric_limits<float>::max())
+        || dval * factor_ >  std::numeric_limits<float>::max())
     {
         GMX_THROW(InvalidInputError("Invalid value: '" + value
                                     + "'; it causes an overflow/underflow"));
index 86c983c955f809758ac25a35d2e40a2ac481ff9e..fc280a5eb48615a347b8a34358ccb818ce7f125f 100644 (file)
@@ -102,7 +102,7 @@ char *fgets2(char *line, int n, FILE *stream)
          * or because of n being too small.
          * Since both cases occur very infrequently, we can check for EOF.
          */
-        if (!gmx_eof(stream))
+        if (!feof(stream))
         {
             gmx_fatal(FARGS, "An input file contains a line longer than %d characters, while the buffer passed to fgets2 has size %d. The line starts with: '%20.20s'", n, n, line);
         }