Merge branch release-2018 into release-2019
[alexxy/gromacs.git] / src / gromacs / mdlib / force.cpp
index c939c4e8a7f13553ec52fde621a8812dc35727b7..5e089e79f65563ddc4611c45fc82f4becd5f49be 100644 (file)
 
 #include "config.h"
 
-#include <assert.h>
-#include <string.h>
-
+#include <cassert>
 #include <cmath>
+#include <cstring>
 
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/listed-forces/listed-forces.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vecdump.h"
+#include "gromacs/mdlib/force_flags.h"
 #include "gromacs/mdlib/forcerec-threading.h"
-#include "gromacs/mdlib/genborn.h"
 #include "gromacs/mdlib/mdrun.h"
 #include "gromacs/mdlib/ns.h"
 #include "gromacs/mdlib/qmmm.h"
+#include "gromacs/mdlib/rf_util.h"
+#include "gromacs/mdlib/wall.h"
 #include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/enerdata.h"
 #include "gromacs/mdtypes/forceoutput.h"
+#include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
-void ns(FILE              *fp,
-        t_forcerec        *fr,
-        matrix             box,
-        gmx_groups_t      *groups,
-        gmx_localtop_t    *top,
-        t_mdatoms         *md,
-        t_commrec         *cr,
-        t_nrnb            *nrnb,
-        gmx_bool           bFillGrid)
+void ns(FILE               *fp,
+        t_forcerec         *fr,
+        matrix              box,
+        const gmx_groups_t *groups,
+        gmx_localtop_t     *top,
+        const t_mdatoms    *md,
+        const t_commrec    *cr,
+        t_nrnb             *nrnb,
+        gmx_bool            bFillGrid)
 {
     int     nsearch;
 
@@ -135,26 +138,28 @@ static void reduceEwaldThreadOuput(int nthreads, ewald_corr_thread_t *ewc_t)
     }
 }
 
-void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
-                       t_idef     *idef,    t_commrec  *cr,
-                       t_nrnb     *nrnb,    gmx_wallcycle_t wcycle,
-                       t_mdatoms  *md,
-                       rvec       x[],      history_t  *hist,
-                       rvec      *forceForUseWithShiftForces,
+void do_force_lowlevel(t_forcerec           *fr,
+                       const t_inputrec     *ir,
+                       const t_idef         *idef,
+                       const t_commrec      *cr,
+                       const gmx_multisim_t *ms,
+                       t_nrnb               *nrnb,
+                       gmx_wallcycle_t       wcycle,
+                       const t_mdatoms      *md,
+                       rvec                  x[],
+                       history_t            *hist,
+                       rvec                 *forceForUseWithShiftForces,
                        gmx::ForceWithVirial *forceWithVirial,
-                       gmx_enerdata_t *enerd,
-                       t_fcdata   *fcd,
-                       gmx_localtop_t *top,
-                       gmx_genborn_t *born,
-                       gmx_bool       bBornRadii,
-                       matrix     box,
-                       t_lambda   *fepvals,
-                       real       *lambda,
-                       t_graph    *graph,
-                       t_blocka   *excl,
-                       rvec       mu_tot[],
-                       int        flags,
-                       float     *cycles_pme)
+                       gmx_enerdata_t       *enerd,
+                       t_fcdata             *fcd,
+                       matrix                box,
+                       t_lambda             *fepvals,
+                       real                 *lambda,
+                       const t_graph        *graph,
+                       const t_blocka       *excl,
+                       rvec                  mu_tot[],
+                       int                   flags,
+                       float                *cycles_pme)
 {
     int         i, j;
     int         donb_flags;
@@ -196,30 +201,12 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
     if (ir->nwall)
     {
         /* foreign lambda component for walls */
-        real dvdl_walls = do_walls(ir, fr, box, md, x, forceForUseWithShiftForces, lambda[efptVDW],
+        real dvdl_walls = do_walls(*ir, *fr, box, *md, x,
+                                   forceWithVirial, lambda[efptVDW],
                                    enerd->grpp.ener[egLJSR], nrnb);
         enerd->dvdl_lin[efptVDW] += dvdl_walls;
     }
 
-    /* If doing GB, reset dvda and calculate the Born radii */
-    if (ir->implicit_solvent)
-    {
-        wallcycle_sub_start(wcycle, ewcsNONBONDED);
-
-        for (i = 0; i < born->nr; i++)
-        {
-            fr->dvda[i] = 0;
-        }
-
-        if (bBornRadii)
-        {
-            calc_gb_rad(cr, fr, ir, top, x, fr->gblist, born, md, nrnb);
-        }
-
-        wallcycle_sub_stop(wcycle, ewcsNONBONDED);
-    }
-
-    where();
     /* We only do non-bonded calculation with group scheme here, the verlet
      * calls are done from do_force_cutsVERLET(). */
     if (fr->cutoff_scheme == ecutsGROUP && (flags & GMX_FORCE_NONBONDED))
@@ -270,18 +257,6 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
             }
         }
         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
-        where();
-    }
-
-    /* If we are doing GB, calculate bonded forces and apply corrections
-     * to the solvation forces */
-    /* MRS: Eventually, many need to include free energy contribution here! */
-    if (ir->implicit_solvent)
-    {
-        wallcycle_sub_start(wcycle, ewcsLISTED);
-        calc_gb_forces(cr, md, born, top, x, forceForUseWithShiftForces, fr, idef,
-                       ir->gb_algorithm, ir->sa_algorithm, nrnb, &pbc, graph, enerd);
-        wallcycle_sub_stop(wcycle, ewcsLISTED);
     }
 
 #if GMX_MPI
@@ -359,14 +334,13 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                    TRUE, box);
     }
 
-    do_force_listed(wcycle, box, ir->fepvals, cr,
-                    idef, (const rvec *) x, hist,
+    do_force_listed(wcycle, box, ir->fepvals, cr, ms,
+                    idef, x, hist,
                     forceForUseWithShiftForces, forceWithVirial,
                     fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
-                    DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr,
+                    DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr,
                     flags);
 
-    where();
 
     *cycles_pme = 0;
 
@@ -429,7 +403,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                                            md->sqrt_c6A, md->sqrt_c6B,
                                            md->sigmaA, md->sigmaB,
                                            md->sigma3A, md->sigma3B,
-                                           md->nChargePerturbed || md->nTypePerturbed,
+                                           (md->nChargePerturbed != 0) || (md->nTypePerturbed != 0),
                                            ir->cutoff_scheme != ecutsVERLET,
                                            excl, x, box, mu_tot,
                                            ir->ewald_geometry,
@@ -580,13 +554,13 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
         {
             real dvdl_rf_excl      = 0;
             enerd->term[F_RF_EXCL] =
-                RF_excl_correction(fr, graph, md, excl, x, forceForUseWithShiftForces,
+                RF_excl_correction(fr, graph, md, excl, DOMAINDECOMP(cr),
+                                   x, forceForUseWithShiftForces,
                                    fr->fshift, &pbc, lambda[efptCOUL], &dvdl_rf_excl);
 
             enerd->dvdl_lin[efptCOUL] += dvdl_rf_excl;
         }
     }
-    where();
 
     if (debug)
     {
@@ -680,7 +654,7 @@ void destroy_enerdata(gmx_enerdata_t *enerd)
     }
 }
 
-static real sum_v(int n, real v[])
+static real sum_v(int n, const real v[])
 {
     real t;
     int  i;
@@ -703,8 +677,6 @@ void sum_epot(gmx_grppairener_t *grpp, real *epot)
     epot[F_LJ]       = sum_v(grpp->nener, grpp->ener[egLJSR]);
     epot[F_LJ14]     = sum_v(grpp->nener, grpp->ener[egLJ14]);
     epot[F_COUL14]   = sum_v(grpp->nener, grpp->ener[egCOUL14]);
-    /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
-    epot[F_GBPOL]   += sum_v(grpp->nener, grpp->ener[egGB]);
 
 /* lattice part of LR doesnt belong to any group
  * and has been added earlier
@@ -793,7 +765,7 @@ void sum_dhdl(gmx_enerdata_t *enerd, gmx::ArrayRef<const real> lambda, t_lambda
 
         double &enerpart_lambda = enerd->enerpart_lambda[i + 1];
 
-        for (size_t j = 0; j < lambda.size(); j++)
+        for (gmx::index j = 0; j < lambda.size(); j++)
         {
             /* Note that this loop is over all dhdl components, not just the separated ones */
             const double dlam  = fepvals->all_lambda[j][i] - lambda[j];