Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / mdlib / force.c
index d4e6445491e36fe2f10bff1cf360d7b6d8ae8050..4313d3e75ef729035af5e0600472183e92fadeb8 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "gromacs/legacyheaders/force.h"
+
+#include "config.h"
 
+#include <assert.h>
 #include <math.h>
 #include <string.h>
-#include <assert.h>
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "macros.h"
-#include "gromacs/utility/smalloc.h"
-#include "macros.h"
-#include "physics.h"
-#include "force.h"
-#include "nonbonded.h"
-#include "names.h"
-#include "network.h"
-#include "pbc.h"
-#include "ns.h"
-#include "nrnb.h"
-#include "bondf.h"
-#include "mshift.h"
-#include "txtdump.h"
-#include "coulomb.h"
-#include "pme.h"
-#include "mdrun.h"
-#include "domdec.h"
-#include "qmmm.h"
-#include "gmx_omp_nthreads.h"
 
+#include "gromacs/legacyheaders/coulomb.h"
+#include "gromacs/legacyheaders/domdec.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/network.h"
+#include "gromacs/legacyheaders/nonbonded.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/ns.h"
+#include "gromacs/legacyheaders/pme.h"
+#include "gromacs/legacyheaders/qmmm.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/listed-forces/bonded.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
 #include "gromacs/timing/wallcycle.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
 void ns(FILE              *fp,
         t_forcerec        *fr,
@@ -140,13 +141,7 @@ static void reduce_thread_forces(int n, rvec *f,
     }
 }
 
-void gmx_print_sepdvdl(FILE *fplog, const char *s, real v, real dvdlambda)
-{
-    fprintf(fplog, "  %-30s V %12.5e  dVdl %12.5e\n", s, v, dvdlambda);
-}
-
-void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
-                       t_forcerec *fr,      t_inputrec *ir,
+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,
@@ -157,7 +152,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                        t_fcdata   *fcd,
                        gmx_localtop_t *top,
                        gmx_genborn_t *born,
-                       t_atomtypes *atype,
                        gmx_bool       bBornRadii,
                        matrix     box,
                        t_lambda   *fepvals,
@@ -170,7 +164,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
 {
     int         i, j;
     int         donb_flags;
-    gmx_bool    bDoEpot, bSepDVDL, bSB;
+    gmx_bool    bDoEpot, bSB;
     int         pme_flags;
     matrix      boxs;
     rvec        box_size;
@@ -184,8 +178,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
     double  t0 = 0.0, t1, t2, t3; /* time measurement for coarse load balancing */
 #endif
 
-#define PRINT_SEPDVDL(s, v, dvdlambda) if (bSepDVDL) { gmx_print_sepdvdl(fplog, s, v, dvdlambda); }
-
     set_pbc(&pbc, fr->ePBC, box);
 
     /* reset free energy components */
@@ -201,7 +193,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         box_size[i] = box[i][i];
     }
 
-    bSepDVDL = (fr->bSepDVDL && do_per_step(step, ir->nstlog));
     debug_gmx();
 
     /* do QMMM first if requested */
@@ -210,12 +201,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         enerd->term[F_EQM] = calculate_QMMM(cr, x, f, fr);
     }
 
-    if (bSepDVDL)
-    {
-        fprintf(fplog, "Step %s: non-bonded V and dVdl for rank %d:\n",
-                gmx_step_str(step, buf), cr->nodeid);
-    }
-
     /* Call the short range functions all in one go. */
 
 #ifdef GMX_MPI
@@ -233,7 +218,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         /* foreign lambda component for walls */
         real dvdl_walls = do_walls(ir, fr, box, md, x, f, lambda[efptVDW],
                                    enerd->grpp.ener[egLJSR], nrnb);
-        PRINT_SEPDVDL("Walls", 0.0, dvdl_walls);
         enerd->dvdl_lin[efptVDW] += dvdl_walls;
     }
 
@@ -316,10 +300,10 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
     /* MRS: Eventually, many need to include free energy contribution here! */
     if (ir->implicit_solvent)
     {
-        wallcycle_sub_start(wcycle, ewcsBONDED);
+        wallcycle_sub_start(wcycle, ewcsLISTED);
         calc_gb_forces(cr, md, born, top, x, f, fr, idef,
                        ir->gb_algorithm, ir->sa_algorithm, nrnb, &pbc, graph, enerd);
-        wallcycle_sub_stop(wcycle, ewcsBONDED);
+        wallcycle_sub_stop(wcycle, ewcsLISTED);
     }
 
 #ifdef GMX_MPI
@@ -351,24 +335,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         enerd->dvdl_lin[efptCOUL] += dvdl_nb[efptCOUL];
     }
 
-    if (bSepDVDL)
-    {
-        real V_short_range    = 0;
-        real dvdl_short_range = 0;
-
-        for (i = 0; i < enerd->grpp.nener; i++)
-        {
-            V_short_range +=
-                (fr->bBHAM ?
-                 enerd->grpp.ener[egBHAMSR][i] :
-                 enerd->grpp.ener[egLJSR][i])
-                + enerd->grpp.ener[egCOULSR][i] + enerd->grpp.ener[egGB][i];
-        }
-        dvdl_short_range = dvdl_nb[efptVDW] + dvdl_nb[efptCOUL];
-        PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",
-                      V_short_range,
-                      dvdl_short_range);
-    }
     debug_gmx();
 
 
@@ -377,9 +343,9 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         pr_rvecs(debug, 0, "fshift after SR", fr->fshift, SHIFTS);
     }
 
-    /* Shift the coordinates. Must be done before bonded forces and PPPM,
+    /* Shift the coordinates. Must be done before listed forces and PPPM,
      * but is also necessary for SHAKE and update, therefore it can NOT
-     * go when no bonded forces have to be evaluated.
+     * go when no listed forces have to be evaluated.
      */
 
     /* Here sometimes we would not need to shift with NBFonly,
@@ -397,9 +363,9 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
             inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
         }
     }
-    /* Check whether we need to do bondeds or correct for exclusions */
+    /* Check whether we need to do listed interactions or correct for exclusions */
     if (fr->bMolPBC &&
-        ((flags & GMX_FORCE_BONDED)
+        ((flags & GMX_FORCE_LISTED)
          || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype)))
     {
         /* Since all atoms are in the rectangular or triclinic unit-cell,
@@ -409,14 +375,13 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
     }
     debug_gmx();
 
-    if (flags & GMX_FORCE_BONDED)
+    if (flags & GMX_FORCE_LISTED)
     {
-        wallcycle_sub_start(wcycle, ewcsBONDED);
-        calc_bonds(fplog, cr->ms,
-                   idef, x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
-                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,
-                   flags,
-                   fr->bSepDVDL && do_per_step(step, ir->nstlog), step);
+        wallcycle_sub_start(wcycle, ewcsLISTED);
+        calc_bonds(cr->ms,
+                   idef, (const rvec *) x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
+                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,
+                   flags);
 
         /* Check if we have to determine energy differences
          * at foreign lambda's.
@@ -435,7 +400,7 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                 {
                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
                 }
-                calc_bonds_lambda(fplog, idef, x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
+                calc_bonds_lambda(idef, (const rvec *) x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
                                   fcd, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
@@ -443,17 +408,25 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
         }
         debug_gmx();
 
-        wallcycle_sub_stop(wcycle, ewcsBONDED);
+        wallcycle_sub_stop(wcycle, ewcsLISTED);
     }
 
     where();
 
     *cycles_pme = 0;
+    clear_mat(fr->vir_el_recip);
+    clear_mat(fr->vir_lj_recip);
+
+    /* Do long-range electrostatics and/or LJ-PME, including related short-range
+     * corrections.
+     */
     if (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype))
     {
-        real Vlr             = 0, Vcorr = 0;
-        real dvdl_long_range = 0;
-        int  status          = 0;
+        real Vlr               = 0, Vcorr = 0;
+        real dvdl_long_range   = 0;
+        int  status            = 0;
+        real Vlr_q             = 0, Vlr_lj = 0, Vcorr_q = 0, Vcorr_lj = 0;
+        real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;
 
         bSB = (ir->nwall == 2);
         if (bSB)
@@ -462,20 +435,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
             svmul(ir->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
             box_size[ZZ] *= ir->wall_ewald_zfac;
         }
-    }
-
-    /* Do long-range electrostatics and/or LJ-PME, including related short-range
-     * corrections.
-     */
-
-    clear_mat(fr->vir_el_recip);
-    clear_mat(fr->vir_lj_recip);
-
-    if (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype))
-    {
-        real Vlr_q             = 0, Vlr_lj = 0, Vcorr_q = 0, Vcorr_lj = 0;
-        real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;
-        int  status            = 0;
 
         if (EEL_PME_EWALD(fr->eeltype) || EVDW_PME(fr->vdwtype))
         {
@@ -577,15 +536,10 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                                                    fr->vir_el_recip);
             }
 
-            PRINT_SEPDVDL("Ewald excl./charge/dip. corr.", Vcorr_q, dvdl_long_range_correction_q);
-            PRINT_SEPDVDL("Ewald excl. corr. LJ", Vcorr_lj, dvdl_long_range_correction_lj);
             enerd->dvdl_lin[efptCOUL] += dvdl_long_range_correction_q;
             enerd->dvdl_lin[efptVDW]  += dvdl_long_range_correction_lj;
-        }
 
-        if ((EEL_PME(fr->eeltype) || EVDW_PME(fr->vdwtype)))
-        {
-            if (cr->duty & DUTY_PME)
+            if ((EEL_PME(fr->eeltype) || EVDW_PME(fr->vdwtype)) && (cr->duty & DUTY_PME))
             {
                 /* Do reciprocal PME for Coulomb and/or LJ. */
                 assert(fr->n_tpi >= 0);
@@ -657,7 +611,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                                         md->chargeA + md->homenr - fr->n_tpi,
                                         &Vlr_q);
                 }
-                PRINT_SEPDVDL("PME mesh", Vlr_q + Vlr_lj, dvdl_long_range_q+dvdl_long_range_lj);
             }
         }
 
@@ -668,7 +621,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                              box_size, cr, md->homenr,
                              fr->vir_el_recip, fr->ewaldcoeff_q,
                              lambda[efptCOUL], &dvdl_long_range_q, fr->ewald_table);
-            PRINT_SEPDVDL("Ewald long-range", Vlr_q, dvdl_long_range_q);
         }
 
         /* Note that with separate PME nodes we get the real energies later */
@@ -703,8 +655,6 @@ void do_force_lowlevel(FILE       *fplog,   gmx_int64_t step,
                                        fr->fshift, &pbc, lambda[efptCOUL], &dvdl_rf_excl);
 
                 enerd->dvdl_lin[efptCOUL] += dvdl_rf_excl;
-                PRINT_SEPDVDL("RF exclusion correction",
-                              enerd->term[F_RF_EXCL], dvdl_rf_excl);
             }
         }
     }