Avoid using function calls in OpenMP directives

[alexxy/gromacs.git] / src / gromacs / mdlib / force.c

diff --git a/src/gromacs/mdlib/force.c b/src/gromacs/mdlib/force.c
index c328a7eddb81a0a3beceb61b48895704758e517a..d4e6445491e36fe2f10bff1cf360d7b6d8ae8050 100644 (file)
--- a/src/gromacs/mdlib/force.c
+++ b/src/gromacs/mdlib/force.c
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014, 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.
@@ -44,7 +44,7 @@
 #include "sysstuff.h"
 #include "typedefs.h"
 #include "macros.h"
-#include "smalloc.h"
+#include "gromacs/utility/smalloc.h"
 #include "macros.h"
 #include "physics.h"
 #include "force.h"
@@ -61,11 +61,11 @@
 #include "pme.h"
 #include "mdrun.h"
 #include "domdec.h"
-#include "partdec.h"
 #include "qmmm.h"
 #include "gmx_omp_nthreads.h"
 
 #include "gromacs/timing/wallcycle.h"
+#include "gmx_fatal.h"
 
 void ns(FILE              *fp,
         t_forcerec        *fr,
@@ -117,9 +117,11 @@ static void reduce_thread_forces(int n, rvec *f,
                                  int nthreads, f_thread_t *f_t)
 {
     int t, i;
+    int nthreads_loop gmx_unused;
 
     /* This reduction can run over any number of threads */
-#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
+    nthreads_loop = gmx_omp_nthreads_get(emntBonded);
+#pragma omp parallel for num_threads(nthreads_loop) private(t) schedule(static)
     for (i = 0; i < n; i++)
     {
         for (t = 1; t < nthreads; t++)
@@ -210,7 +212,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
 
     if (bSepDVDL)
     {
-        fprintf(fplog, "Step %s: non-bonded V and dVdl for node %d:\n",
+        fprintf(fplog, "Step %s: non-bonded V and dVdl for rank %d:\n",
                 gmx_step_str(step, buf), cr->nodeid);
     }
 
@@ -262,10 +264,15 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         /* Add short-range interactions */
         donb_flags |= GMX_NONBONDED_DO_SR;
 
+        /* Currently all group scheme kernels always calculate (shift-)forces */
         if (flags & GMX_FORCE_FORCES)
         {
             donb_flags |= GMX_NONBONDED_DO_FORCE;
         }
+        if (flags & GMX_FORCE_VIRIAL)
+        {
+            donb_flags |= GMX_NONBONDED_DO_SHIFTFORCE;
+        }
         if (flags & GMX_FORCE_ENERGY)
         {
             donb_flags |= GMX_NONBONDED_DO_POTENTIAL;
@@ -470,7 +477,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;
         int  status            = 0;
 
-        if (EEL_EWALD(fr->eeltype) || EVDW_PME(fr->vdwtype))
+        if (EEL_PME_EWALD(fr->eeltype) || EVDW_PME(fr->vdwtype))
         {
             real dvdl_long_range_correction_q   = 0;
             real dvdl_long_range_correction_lj  = 0;
@@ -531,16 +538,17 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                     }
                     *dvdlt_q  = 0;
                     *dvdlt_lj = 0;
+
                     ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
                                        cr, t, fr,
                                        md->chargeA,
                                        md->nChargePerturbed ? md->chargeB : NULL,
-                                       md->c6A,
-                                       md->nChargePerturbed ? md->c6B : NULL,
+                                       md->sqrt_c6A,
+                                       md->nTypePerturbed ? md->sqrt_c6B : NULL,
                                        md->sigmaA,
-                                       md->nChargePerturbed ? md->sigmaB : NULL,
+                                       md->nTypePerturbed ? md->sigmaB : NULL,
                                        md->sigma3A,
-                                       md->nChargePerturbed ? md->sigma3B : NULL,
+                                       md->nTypePerturbed ? md->sigma3B : NULL,
                                        ir->cutoff_scheme != ecutsVERLET,
                                        excl, x, bSB ? boxs : box, mu_tot,
                                        ir->ewald_geometry,
@@ -562,7 +570,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                 wallcycle_sub_stop(wcycle, ewcsEWALD_CORRECTION);
             }
 
-            if (fr->n_tpi == 0)
+            if (EEL_PME_EWALD(fr->eeltype) && fr->n_tpi == 0)
             {
                 Vcorr_q += ewald_charge_correction(cr, fr, lambda[efptCOUL], box,
                                                    &dvdl_long_range_correction_q,
@@ -583,7 +591,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                 assert(fr->n_tpi >= 0);
                 if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
                 {
-                    pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
+                    pme_flags = GMX_PME_SPREAD | GMX_PME_SOLVE;
                     if (EEL_PME(fr->eeltype))
                     {
                         pme_flags     |= GMX_PME_DO_COULOMB;
@@ -591,11 +599,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                     if (EVDW_PME(fr->vdwtype))
                     {
                         pme_flags |= GMX_PME_DO_LJ;
-                        if (fr->ljpme_combination_rule == eljpmeLB)
-                        {
-                            /*Lorentz-Berthelot Comb. Rules in LJ-PME*/
-                            pme_flags |= GMX_PME_LJ_LB;
-                        }
                     }
                     if (flags & GMX_FORCE_FORCES)
                     {
@@ -612,10 +615,10 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                     }
                     wallcycle_start(wcycle, ewcPMEMESH);
                     status = gmx_pme_do(fr->pmedata,
-                                        md->start, md->homenr - fr->n_tpi,
+                                        0, md->homenr - fr->n_tpi,
                                         x, fr->f_novirsum,
                                         md->chargeA, md->chargeB,
-                                        md->c6A, md->c6B,
+                                        md->sqrt_c6A, md->sqrt_c6B,
                                         md->sigmaA, md->sigmaB,
                                         bSB ? boxs : box, cr,
                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
@@ -658,7 +661,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
             }
         }
 
-        if (!EEL_PME(fr->eeltype) && EEL_EWALD(fr->eeltype))
+        if (!EEL_PME(fr->eeltype) && EEL_PME_EWALD(fr->eeltype))
         {
             Vlr_q = do_ewald(ir, x, fr->f_novirsum,
                              md->chargeA, md->chargeB,
@@ -667,11 +670,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                              lambda[efptCOUL], &dvdl_long_range_q, fr->ewald_table);
             PRINT_SEPDVDL("Ewald long-range", Vlr_q, dvdl_long_range_q);
         }
-        else if (!EEL_EWALD(fr->eeltype))
-        {
-            gmx_fatal(FARGS, "No such electrostatics method implemented %s",
-                      eel_names[fr->eeltype]);
-        }
+
         /* Note that with separate PME nodes we get the real energies later */
         enerd->dvdl_lin[efptCOUL] += dvdl_long_range_q;
         enerd->dvdl_lin[efptVDW]  += dvdl_long_range_lj;
@@ -727,7 +726,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         fr->t_wait += t3-t2;
         if (fr->timesteps == 11)
         {
-            fprintf(stderr, "* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
+            fprintf(stderr, "* PP load balancing info: rank %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
                     cr->nodeid, gmx_step_str(fr->timesteps, buf),
                     100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
                     (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);