src/gromacs/mdlib/force.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include <assert.h>
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "physics.h"
  49 #include "force.h"
  50 #include "nonbonded.h"
  51 #include "names.h"
  52 #include "network.h"
  53 #include "pbc.h"
  54 #include "ns.h"
  55 #include "nrnb.h"
  56 #include "bondf.h"
  57 #include "mshift.h"
  58 #include "txtdump.h"
  59 #include "coulomb.h"
  60 #include "pme.h"
  61 #include "mdrun.h"
  62 #include "domdec.h"
  63 #include "partdec.h"
  64 #include "qmmm.h"
  65 #include "gmx_omp_nthreads.h"
  66
  67
  68 void ns(FILE              *fp,
  69         t_forcerec        *fr,
  70         matrix             box,
  71         gmx_groups_t      *groups,
  72         gmx_localtop_t    *top,
  73         t_mdatoms         *md,
  74         t_commrec         *cr,
  75         t_nrnb            *nrnb,
  76         gmx_bool           bFillGrid,
  77         gmx_bool           bDoLongRangeNS)
  78 {
  79     char   *ptr;
  80     int     nsearch;
  81
  82
  83     if (!fr->ns.nblist_initialized)
  84     {
  85         init_neighbor_list(fp, fr, md->homenr);
  86     }
  87
  88     if (fr->bTwinRange)
  89     {
  90         fr->nlr = 0;
  91     }
  92
  93     nsearch = search_neighbours(fp, fr, box, top, groups, cr, nrnb, md,
  94                                 bFillGrid, bDoLongRangeNS);
  95     if (debug)
  96     {
  97         fprintf(debug, "nsearch = %d\n", nsearch);
  98     }
  99
 100     /* Check whether we have to do dynamic load balancing */
 101     /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
 102        count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
 103        &(top->idef),opts->ngener);
 104      */
 105     if (fr->ns.dump_nl > 0)
 106     {
 107         dump_nblist(fp, cr, fr, fr->ns.dump_nl);
 108     }
 109 }
 110
 111 static void reduce_thread_forces(int n, rvec *f,
 112                                  tensor vir,
 113                                  real *Vcorr,
 114                                  int efpt_ind, real *dvdl,
 115                                  int nthreads, f_thread_t *f_t)
 116 {
 117     int t, i;
 118
 119     /* This reduction can run over any number of threads */
 120 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
 121     for (i = 0; i < n; i++)
 122     {
 123         for (t = 1; t < nthreads; t++)
 124         {
 125             rvec_inc(f[i], f_t[t].f[i]);
 126         }
 127     }
 128     for (t = 1; t < nthreads; t++)
 129     {
 130         *Vcorr += f_t[t].Vcorr;
 131         *dvdl  += f_t[t].dvdl[efpt_ind];
 132         m_add(vir, f_t[t].vir, vir);
 133     }
 134 }
 135
 136 void gmx_print_sepdvdl(FILE *fplog, const char *s, real v, real dvdlambda)
 137 {
 138     fprintf(fplog, "  %-30s V %12.5e  dVdl %12.5e\n", s, v, dvdlambda);
 139 }
 140
 141 void do_force_lowlevel(FILE       *fplog,   gmx_large_int_t step,
 142                        t_forcerec *fr,      t_inputrec *ir,
 143                        t_idef     *idef,    t_commrec  *cr,
 144                        t_nrnb     *nrnb,    gmx_wallcycle_t wcycle,
 145                        t_mdatoms  *md,
 146                        rvec       x[],      history_t  *hist,
 147                        rvec       f[],
 148                        rvec       f_longrange[],
 149                        gmx_enerdata_t *enerd,
 150                        t_fcdata   *fcd,
 151                        gmx_localtop_t *top,
 152                        gmx_genborn_t *born,
 153                        t_atomtypes *atype,
 154                        gmx_bool       bBornRadii,
 155                        matrix     box,
 156                        t_lambda   *fepvals,
 157                        real       *lambda,
 158                        t_graph    *graph,
 159                        t_blocka   *excl,
 160                        rvec       mu_tot[],
 161                        int        flags,
 162                        float      *cycles_pme)
 163 {
 164     int         i, j, status;
 165     int         donb_flags;
 166     gmx_bool    bDoEpot, bSepDVDL, bSB;
 167     int         pme_flags;
 168     matrix      boxs;
 169     rvec        box_size;
 170     real        Vsr, Vlr, Vcorr = 0;
 171     t_pbc       pbc;
 172     real        dvdgb;
 173     char        buf[22];
 174     double      clam_i, vlam_i;
 175     real        dvdl_dum[efptNR], dvdl, dvdl_nb[efptNR], lam_i[efptNR];
 176     real        dvdlsum;
 177
 178 #ifdef GMX_MPI
 179     double  t0 = 0.0, t1, t2, t3; /* time measurement for coarse load balancing */
 180 #endif
 181
 182 #define PRINT_SEPDVDL(s, v, dvdlambda) if (bSepDVDL) { gmx_print_sepdvdl(fplog, s, v, dvdlambda); }
 183
 184     set_pbc(&pbc, fr->ePBC, box);
 185
 186     /* reset free energy components */
 187     for (i = 0; i < efptNR; i++)
 188     {
 189         dvdl_nb[i]  = 0;
 190         dvdl_dum[i] = 0;
 191     }
 192
 193     /* Reset box */
 194     for (i = 0; (i < DIM); i++)
 195     {
 196         box_size[i] = box[i][i];
 197     }
 198
 199     bSepDVDL = (fr->bSepDVDL && do_per_step(step, ir->nstlog));
 200     debug_gmx();
 201
 202     /* do QMMM first if requested */
 203     if (fr->bQMMM)
 204     {
 205         enerd->term[F_EQM] = calculate_QMMM(cr, x, f, fr);
 206     }
 207
 208     if (bSepDVDL)
 209     {
 210         fprintf(fplog, "Step %s: non-bonded V and dVdl for node %d:\n",
 211                 gmx_step_str(step, buf), cr->nodeid);
 212     }
 213
 214     /* Call the short range functions all in one go. */
 215
 216 #ifdef GMX_MPI
 217     /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
 218 #define TAKETIME FALSE
 219     if (TAKETIME)
 220     {
 221         MPI_Barrier(cr->mpi_comm_mygroup);
 222         t0 = MPI_Wtime();
 223     }
 224 #endif
 225
 226     if (ir->nwall)
 227     {
 228         /* foreign lambda component for walls */
 229         dvdl = do_walls(ir, fr, box, md, x, f, lambda[efptVDW],
 230                         enerd->grpp.ener[egLJSR], nrnb);
 231         PRINT_SEPDVDL("Walls", 0.0, dvdl);
 232         enerd->dvdl_lin[efptVDW] += dvdl;
 233     }
 234
 235     /* If doing GB, reset dvda and calculate the Born radii */
 236     if (ir->implicit_solvent)
 237     {
 238         wallcycle_sub_start(wcycle, ewcsNONBONDED);
 239
 240         for (i = 0; i < born->nr; i++)
 241         {
 242             fr->dvda[i] = 0;
 243         }
 244
 245         if (bBornRadii)
 246         {
 247             calc_gb_rad(cr, fr, ir, top, x, &(fr->gblist), born, md, nrnb);
 248         }
 249
 250         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
 251     }
 252
 253     where();
 254     /* We only do non-bonded calculation with group scheme here, the verlet
 255      * calls are done from do_force_cutsVERLET(). */
 256     if (fr->cutoff_scheme == ecutsGROUP && (flags & GMX_FORCE_NONBONDED))
 257     {
 258         donb_flags = 0;
 259         /* Add short-range interactions */
 260         donb_flags |= GMX_NONBONDED_DO_SR;
 261
 262         if (flags & GMX_FORCE_FORCES)
 263         {
 264             donb_flags |= GMX_NONBONDED_DO_FORCE;
 265         }
 266         if (flags & GMX_FORCE_ENERGY)
 267         {
 268             donb_flags |= GMX_NONBONDED_DO_POTENTIAL;
 269         }
 270         if (flags & GMX_FORCE_DO_LR)
 271         {
 272             donb_flags |= GMX_NONBONDED_DO_LR;
 273         }
 274
 275         wallcycle_sub_start(wcycle, ewcsNONBONDED);
 276         do_nonbonded(fr, x, f, f_longrange, md, excl,
 277                      &enerd->grpp, nrnb,
 278                      lambda, dvdl_nb, -1, -1, donb_flags);
 279
 280         /* If we do foreign lambda and we have soft-core interactions
 281          * we have to recalculate the (non-linear) energies contributions.
 282          */
 283         if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL) && fepvals->sc_alpha != 0)
 284         {
 285             for (i = 0; i < enerd->n_lambda; i++)
 286             {
 287                 for (j = 0; j < efptNR; j++)
 288                 {
 289                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 290                 }
 291                 reset_foreign_enerdata(enerd);
 292                 do_nonbonded(fr, x, f, f_longrange, md, excl,
 293                              &(enerd->foreign_grpp), nrnb,
 294                              lam_i, dvdl_dum, -1, -1,
 295                              (donb_flags & ~GMX_NONBONDED_DO_FORCE) | GMX_NONBONDED_DO_FOREIGNLAMBDA);
 296                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 297                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 298             }
 299         }
 300         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
 301         where();
 302     }
 303
 304     /* If we are doing GB, calculate bonded forces and apply corrections
 305      * to the solvation forces */
 306     /* MRS: Eventually, many need to include free energy contribution here! */
 307     if (ir->implicit_solvent)
 308     {
 309         wallcycle_sub_start(wcycle, ewcsBONDED);
 310         calc_gb_forces(cr, md, born, top, x, f, fr, idef,
 311                        ir->gb_algorithm, ir->sa_algorithm, nrnb, &pbc, graph, enerd);
 312         wallcycle_sub_stop(wcycle, ewcsBONDED);
 313     }
 314
 315 #ifdef GMX_MPI
 316     if (TAKETIME)
 317     {
 318         t1          = MPI_Wtime();
 319         fr->t_fnbf += t1-t0;
 320     }
 321 #endif
 322
 323     if (fepvals->sc_alpha != 0)
 324     {
 325         enerd->dvdl_nonlin[efptVDW] += dvdl_nb[efptVDW];
 326     }
 327     else
 328     {
 329         enerd->dvdl_lin[efptVDW] += dvdl_nb[efptVDW];
 330     }
 331
 332     if (fepvals->sc_alpha != 0)
 333
 334     /* even though coulomb part is linear, we already added it, beacuse we
 335        need to go through the vdw calculation anyway */
 336     {
 337         enerd->dvdl_nonlin[efptCOUL] += dvdl_nb[efptCOUL];
 338     }
 339     else
 340     {
 341         enerd->dvdl_lin[efptCOUL] += dvdl_nb[efptCOUL];
 342     }
 343
 344     Vsr = 0;
 345     if (bSepDVDL)
 346     {
 347         for (i = 0; i < enerd->grpp.nener; i++)
 348         {
 349             Vsr +=
 350                 (fr->bBHAM ?
 351                  enerd->grpp.ener[egBHAMSR][i] :
 352                  enerd->grpp.ener[egLJSR][i])
 353                 + enerd->grpp.ener[egCOULSR][i] + enerd->grpp.ener[egGB][i];
 354         }
 355         dvdlsum = dvdl_nb[efptVDW] + dvdl_nb[efptCOUL];
 356         PRINT_SEPDVDL("VdW and Coulomb SR particle-p.", Vsr, dvdlsum);
 357     }
 358     debug_gmx();
 359
 360
 361     if (debug)
 362     {
 363         pr_rvecs(debug, 0, "fshift after SR", fr->fshift, SHIFTS);
 364     }
 365
 366     /* Shift the coordinates. Must be done before bonded forces and PPPM,
 367      * but is also necessary for SHAKE and update, therefore it can NOT
 368      * go when no bonded forces have to be evaluated.
 369      */
 370
 371     /* Here sometimes we would not need to shift with NBFonly,
 372      * but we do so anyhow for consistency of the returned coordinates.
 373      */
 374     if (graph)
 375     {
 376         shift_self(graph, box, x);
 377         if (TRICLINIC(box))
 378         {
 379             inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
 380         }
 381         else
 382         {
 383             inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
 384         }
 385     }
 386     /* Check whether we need to do bondeds or correct for exclusions */
 387     if (fr->bMolPBC &&
 388         ((flags & GMX_FORCE_BONDED)
 389          || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))
 390     {
 391         /* Since all atoms are in the rectangular or triclinic unit-cell,
 392          * only single box vector shifts (2 in x) are required.
 393          */
 394         set_pbc_dd(&pbc, fr->ePBC, cr->dd, TRUE, box);
 395     }
 396     debug_gmx();
 397
 398     if (flags & GMX_FORCE_BONDED)
 399     {
 400         wallcycle_sub_start(wcycle, ewcsBONDED);
 401         calc_bonds(fplog, cr->ms,
 402                    idef, x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
 403                    DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,
 404                    flags,
 405                    fr->bSepDVDL && do_per_step(step, ir->nstlog), step);
 406
 407         /* Check if we have to determine energy differences
 408          * at foreign lambda's.
 409          */
 410         if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL) &&
 411             idef->ilsort != ilsortNO_FE)
 412         {
 413             if (idef->ilsort != ilsortFE_SORTED)
 414             {
 415                 gmx_incons("The bonded interactions are not sorted for free energy");
 416             }
 417             for (i = 0; i < enerd->n_lambda; i++)
 418             {
 419                 reset_foreign_enerdata(enerd);
 420                 for (j = 0; j < efptNR; j++)
 421                 {
 422                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 423                 }
 424                 calc_bonds_lambda(fplog, idef, x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 425                                   fcd, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
 426                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 427                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 428             }
 429         }
 430         debug_gmx();
 431
 432         wallcycle_sub_stop(wcycle, ewcsBONDED);
 433     }
 434
 435     where();
 436
 437     *cycles_pme = 0;
 438     if (EEL_FULL(fr->eeltype))
 439     {
 440         bSB = (ir->nwall == 2);
 441         if (bSB)
 442         {
 443             copy_mat(box, boxs);
 444             svmul(ir->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
 445             box_size[ZZ] *= ir->wall_ewald_zfac;
 446         }
 447
 448         clear_mat(fr->vir_el_recip);
 449
 450         if (fr->bEwald)
 451         {
 452             Vcorr = 0;
 453             dvdl  = 0;
 454
 455             /* With the Verlet scheme exclusion forces are calculated
 456              * in the non-bonded kernel.
 457              */
 458             /* The TPI molecule does not have exclusions with the rest
 459              * of the system and no intra-molecular PME grid contributions
 460              * will be calculated in gmx_pme_calc_energy.
 461              */
 462             if ((ir->cutoff_scheme == ecutsGROUP && fr->n_tpi == 0) ||
 463                 ir->ewald_geometry != eewg3D ||
 464                 ir->epsilon_surface != 0)
 465             {
 466                 int nthreads, t;
 467
 468                 wallcycle_sub_start(wcycle, ewcsEWALD_CORRECTION);
 469
 470                 if (fr->n_tpi > 0)
 471                 {
 472                     gmx_fatal(FARGS, "TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
 473                 }
 474
 475                 nthreads = gmx_omp_nthreads_get(emntBonded);
 476 #pragma omp parallel for num_threads(nthreads) schedule(static)
 477                 for (t = 0; t < nthreads; t++)
 478                 {
 479                     int     s, e, i;
 480                     rvec   *fnv;
 481                     tensor *vir;
 482                     real   *Vcorrt, *dvdlt;
 483                     if (t == 0)
 484                     {
 485                         fnv    = fr->f_novirsum;
 486                         vir    = &fr->vir_el_recip;
 487                         Vcorrt = &Vcorr;
 488                         dvdlt  = &dvdl;
 489                     }
 490                     else
 491                     {
 492                         fnv    = fr->f_t[t].f;
 493                         vir    = &fr->f_t[t].vir;
 494                         Vcorrt = &fr->f_t[t].Vcorr;
 495                         dvdlt  = &fr->f_t[t].dvdl[efptCOUL];
 496                         for (i = 0; i < fr->natoms_force; i++)
 497                         {
 498                             clear_rvec(fnv[i]);
 499                         }
 500                         clear_mat(*vir);
 501                     }
 502                     *dvdlt  = 0;
 503                     *Vcorrt =
 504                         ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
 505                                            cr, t, fr,
 506                                            md->chargeA,
 507                                            md->nChargePerturbed ? md->chargeB : NULL,
 508                                            ir->cutoff_scheme != ecutsVERLET,
 509                                            excl, x, bSB ? boxs : box, mu_tot,
 510                                            ir->ewald_geometry,
 511                                            ir->epsilon_surface,
 512                                            fnv, *vir,
 513                                            lambda[efptCOUL], dvdlt);
 514                 }
 515                 if (nthreads > 1)
 516                 {
 517                     reduce_thread_forces(fr->natoms_force, fr->f_novirsum,
 518                                          fr->vir_el_recip,
 519                                          &Vcorr, efptCOUL, &dvdl,
 520                                          nthreads, fr->f_t);
 521                 }
 522
 523                 wallcycle_sub_stop(wcycle, ewcsEWALD_CORRECTION);
 524             }
 525
 526             if (fr->n_tpi == 0)
 527             {
 528                 Vcorr += ewald_charge_correction(cr, fr, lambda[efptCOUL], box,
 529                                                  &dvdl, fr->vir_el_recip);
 530             }
 531
 532             PRINT_SEPDVDL("Ewald excl./charge/dip. corr.", Vcorr, dvdl);
 533             enerd->dvdl_lin[efptCOUL] += dvdl;
 534         }
 535
 536         status = 0;
 537         Vlr    = 0;
 538         dvdl   = 0;
 539         switch (fr->eeltype)
 540         {
 541             case eelPME:
 542             case eelPMESWITCH:
 543             case eelPMEUSER:
 544             case eelPMEUSERSWITCH:
 545             case eelP3M_AD:
 546                 if (cr->duty & DUTY_PME)
 547                 {
 548                     assert(fr->n_tpi >= 0);
 549                     if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
 550                     {
 551                         pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
 552                         if (flags & GMX_FORCE_FORCES)
 553                         {
 554                             pme_flags |= GMX_PME_CALC_F;
 555                         }
 556                         if (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY))
 557                         {
 558                             pme_flags |= GMX_PME_CALC_ENER_VIR;
 559                         }
 560                         if (fr->n_tpi > 0)
 561                         {
 562                             /* We don't calculate f, but we do want the potential */
 563                             pme_flags |= GMX_PME_CALC_POT;
 564                         }
 565                         wallcycle_start(wcycle, ewcPMEMESH);
 566                         status = gmx_pme_do(fr->pmedata,
 567                                             md->start, md->homenr - fr->n_tpi,
 568                                             x, fr->f_novirsum,
 569                                             md->chargeA, md->chargeB,
 570                                             bSB ? boxs : box, cr,
 571                                             DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
 572                                             DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
 573                                             nrnb, wcycle,
 574                                             fr->vir_el_recip, fr->ewaldcoeff,
 575                                             &Vlr, lambda[efptCOUL], &dvdl,
 576                                             pme_flags);
 577                         *cycles_pme = wallcycle_stop(wcycle, ewcPMEMESH);
 578
 579                         /* We should try to do as little computation after
 580                          * this as possible, because parallel PME synchronizes
 581                          * the nodes, so we want all load imbalance of the rest
 582                          * of the force calculation to be before the PME call.
 583                          * DD load balancing is done on the whole time of
 584                          * the force call (without PME).
 585                          */
 586                     }
 587                     if (fr->n_tpi > 0)
 588                     {
 589                         /* Determine the PME grid energy of the test molecule
 590                          * with the PME grid potential of the other charges.
 591                          */
 592                         gmx_pme_calc_energy(fr->pmedata, fr->n_tpi,
 593                                             x + md->homenr - fr->n_tpi,
 594                                             md->chargeA + md->homenr - fr->n_tpi,
 595                                             &Vlr);
 596                     }
 597                     PRINT_SEPDVDL("PME mesh", Vlr, dvdl);
 598                 }
 599                 break;
 600             case eelEWALD:
 601                 Vlr = do_ewald(ir, x, fr->f_novirsum,
 602                                md->chargeA, md->chargeB,
 603                                box_size, cr, md->homenr,
 604                                fr->vir_el_recip, fr->ewaldcoeff,
 605                                lambda[efptCOUL], &dvdl, fr->ewald_table);
 606                 PRINT_SEPDVDL("Ewald long-range", Vlr, dvdl);
 607                 break;
 608             default:
 609                 gmx_fatal(FARGS, "No such electrostatics method implemented %s",
 610                           eel_names[fr->eeltype]);
 611         }
 612         if (status != 0)
 613         {
 614             gmx_fatal(FARGS, "Error %d in long range electrostatics routine %s",
 615                       status, EELTYPE(fr->eeltype));
 616         }
 617         /* Note that with separate PME nodes we get the real energies later */
 618         enerd->dvdl_lin[efptCOUL] += dvdl;
 619         enerd->term[F_COUL_RECIP]  = Vlr + Vcorr;
 620         if (debug)
 621         {
 622             fprintf(debug, "Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
 623                     Vlr, Vcorr, enerd->term[F_COUL_RECIP]);
 624             pr_rvecs(debug, 0, "vir_el_recip after corr", fr->vir_el_recip, DIM);
 625             pr_rvecs(debug, 0, "fshift after LR Corrections", fr->fshift, SHIFTS);
 626         }
 627     }
 628     else
 629     {
 630         if (EEL_RF(fr->eeltype))
 631         {
 632             /* With the Verlet scheme exclusion forces are calculated
 633              * in the non-bonded kernel.
 634              */
 635             if (ir->cutoff_scheme != ecutsVERLET && fr->eeltype != eelRF_NEC)
 636             {
 637                 dvdl                   = 0;
 638                 enerd->term[F_RF_EXCL] =
 639                     RF_excl_correction(fr, graph, md, excl, x, f,
 640                                        fr->fshift, &pbc, lambda[efptCOUL], &dvdl);
 641             }
 642
 643             enerd->dvdl_lin[efptCOUL] += dvdl;
 644             PRINT_SEPDVDL("RF exclusion correction",
 645                           enerd->term[F_RF_EXCL], dvdl);
 646         }
 647     }
 648     where();
 649     debug_gmx();
 650
 651     if (debug)
 652     {
 653         print_nrnb(debug, nrnb);
 654     }
 655     debug_gmx();
 656
 657 #ifdef GMX_MPI
 658     if (TAKETIME)
 659     {
 660         t2 = MPI_Wtime();
 661         MPI_Barrier(cr->mpi_comm_mygroup);
 662         t3          = MPI_Wtime();
 663         fr->t_wait += t3-t2;
 664         if (fr->timesteps == 11)
 665         {
 666             fprintf(stderr, "* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
 667                     cr->nodeid, gmx_step_str(fr->timesteps, buf),
 668                     100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
 669                     (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);
 670         }
 671         fr->timesteps++;
 672     }
 673 #endif
 674
 675     if (debug)
 676     {
 677         pr_rvecs(debug, 0, "fshift after bondeds", fr->fshift, SHIFTS);
 678     }
 679
 680 }
 681
 682 void init_enerdata(int ngener, int n_lambda, gmx_enerdata_t *enerd)
 683 {
 684     int i, n2;
 685
 686     for (i = 0; i < F_NRE; i++)
 687     {
 688         enerd->term[i]         = 0;
 689         enerd->foreign_term[i] = 0;
 690     }
 691
 692
 693     for (i = 0; i < efptNR; i++)
 694     {
 695         enerd->dvdl_lin[i]     = 0;
 696         enerd->dvdl_nonlin[i]  = 0;
 697     }
 698
 699     n2 = ngener*ngener;
 700     if (debug)
 701     {
 702         fprintf(debug, "Creating %d sized group matrix for energies\n", n2);
 703     }
 704     enerd->grpp.nener         = n2;
 705     enerd->foreign_grpp.nener = n2;
 706     for (i = 0; (i < egNR); i++)
 707     {
 708         snew(enerd->grpp.ener[i], n2);
 709         snew(enerd->foreign_grpp.ener[i], n2);
 710     }
 711
 712     if (n_lambda)
 713     {
 714         enerd->n_lambda = 1 + n_lambda;
 715         snew(enerd->enerpart_lambda, enerd->n_lambda);
 716     }
 717     else
 718     {
 719         enerd->n_lambda = 0;
 720     }
 721 }
 722
 723 void destroy_enerdata(gmx_enerdata_t *enerd)
 724 {
 725     int i;
 726
 727     for (i = 0; (i < egNR); i++)
 728     {
 729         sfree(enerd->grpp.ener[i]);
 730     }
 731
 732     for (i = 0; (i < egNR); i++)
 733     {
 734         sfree(enerd->foreign_grpp.ener[i]);
 735     }
 736
 737     if (enerd->n_lambda)
 738     {
 739         sfree(enerd->enerpart_lambda);
 740     }
 741 }
 742
 743 static real sum_v(int n, real v[])
 744 {
 745     real t;
 746     int  i;
 747
 748     t = 0.0;
 749     for (i = 0; (i < n); i++)
 750     {
 751         t = t + v[i];
 752     }
 753
 754     return t;
 755 }
 756
 757 void sum_epot(gmx_grppairener_t *grpp, real *epot)
 758 {
 759     int i;
 760
 761     /* Accumulate energies */
 762     epot[F_COUL_SR]  = sum_v(grpp->nener, grpp->ener[egCOULSR]);
 763     epot[F_LJ]       = sum_v(grpp->nener, grpp->ener[egLJSR]);
 764     epot[F_LJ14]     = sum_v(grpp->nener, grpp->ener[egLJ14]);
 765     epot[F_COUL14]   = sum_v(grpp->nener, grpp->ener[egCOUL14]);
 766     epot[F_COUL_LR]  = sum_v(grpp->nener, grpp->ener[egCOULLR]);
 767     epot[F_LJ_LR]    = sum_v(grpp->nener, grpp->ener[egLJLR]);
 768     /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
 769     epot[F_GBPOL]   += sum_v(grpp->nener, grpp->ener[egGB]);
 770
 771 /* lattice part of LR doesnt belong to any group
 772  * and has been added earlier
 773  */
 774     epot[F_BHAM]     = sum_v(grpp->nener, grpp->ener[egBHAMSR]);
 775     epot[F_BHAM_LR]  = sum_v(grpp->nener, grpp->ener[egBHAMLR]);
 776
 777     epot[F_EPOT] = 0;
 778     for (i = 0; (i < F_EPOT); i++)
 779     {
 780         if (i != F_DISRESVIOL && i != F_ORIRESDEV)
 781         {
 782             epot[F_EPOT] += epot[i];
 783         }
 784     }
 785 }
 786
 787 void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
 788 {
 789     int    i, j, index;
 790     double dlam;
 791
 792     enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW];  /* include dispersion correction */
 793     enerd->term[F_DVDL]       = 0.0;
 794     for (i = 0; i < efptNR; i++)
 795     {
 796         if (fepvals->separate_dvdl[i])
 797         {
 798             /* could this be done more readably/compactly? */
 799             switch (i)
 800             {
 801                 case (efptMASS):
 802                     index = F_DKDL;
 803                     break;
 804                 case (efptCOUL):
 805                     index = F_DVDL_COUL;
 806                     break;
 807                 case (efptVDW):
 808                     index = F_DVDL_VDW;
 809                     break;
 810                 case (efptBONDED):
 811                     index = F_DVDL_BONDED;
 812                     break;
 813                 case (efptRESTRAINT):
 814                     index = F_DVDL_RESTRAINT;
 815                     break;
 816                 default:
 817                     index = F_DVDL;
 818                     break;
 819             }
 820             enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 821             if (debug)
 822             {
 823                 fprintf(debug, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
 824                         efpt_names[i], i, enerd->term[index], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
 825             }
 826         }
 827         else
 828         {
 829             enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 830             if (debug)
 831             {
 832                 fprintf(debug, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
 833                         efpt_names[0], i, enerd->term[F_DVDL], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
 834             }
 835         }
 836     }
 837
 838     /* Notes on the foreign lambda free energy difference evaluation:
 839      * Adding the potential and ekin terms that depend linearly on lambda
 840      * as delta lam * dvdl to the energy differences is exact.
 841      * For the constraints this is not exact, but we have no other option
 842      * without literally changing the lengths and reevaluating the energies at each step.
 843      * (try to remedy this post 4.6 - MRS)
 844      * For the non-bonded LR term we assume that the soft-core (if present)
 845      * no longer affects the energy beyond the short-range cut-off,
 846      * which is a very good approximation (except for exotic settings).
 847      * (investigate how to overcome this post 4.6 - MRS)
 848      */
 849     if (fepvals->separate_dvdl[efptBONDED])
 850     {
 851         enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
 852     }
 853     else
 854     {
 855         enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
 856     }
 857     enerd->term[F_DVDL_CONSTR] = 0;
 858
 859     for (i = 0; i < fepvals->n_lambda; i++)
 860     {
 861         /* note we are iterating over fepvals here!
 862            For the current lam, dlam = 0 automatically,
 863            so we don't need to add anything to the
 864            enerd->enerpart_lambda[0] */
 865
 866         /* we don't need to worry about dvdl_lin contributions to dE at
 867            current lambda, because the contributions to the current
 868            lambda are automatically zeroed */
 869
 870         for (j = 0; j < efptNR; j++)
 871         {
 872             /* Note that this loop is over all dhdl components, not just the separated ones */
 873             dlam = (fepvals->all_lambda[j][i]-lambda[j]);
 874             enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
 875             if (debug)
 876             {
 877                 fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
 878                         fepvals->all_lambda[j][i], efpt_names[j],
 879                         (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),
 880                         dlam, enerd->dvdl_lin[j]);
 881             }
 882         }
 883     }
 884 }
 885
 886
 887 void reset_foreign_enerdata(gmx_enerdata_t *enerd)
 888 {
 889     int  i, j;
 890
 891     /* First reset all foreign energy components.  Foreign energies always called on
 892        neighbor search steps */
 893     for (i = 0; (i < egNR); i++)
 894     {
 895         for (j = 0; (j < enerd->grpp.nener); j++)
 896         {
 897             enerd->foreign_grpp.ener[i][j] = 0.0;
 898         }
 899     }
 900
 901     /* potential energy components */
 902     for (i = 0; (i <= F_EPOT); i++)
 903     {
 904         enerd->foreign_term[i] = 0.0;
 905     }
 906 }
 907
 908 void reset_enerdata(t_forcerec *fr, gmx_bool bNS,
 909                     gmx_enerdata_t *enerd,
 910                     gmx_bool bMaster)
 911 {
 912     gmx_bool bKeepLR;
 913     int      i, j;
 914
 915     /* First reset all energy components, except for the long range terms
 916      * on the master at non neighbor search steps, since the long range
 917      * terms have already been summed at the last neighbor search step.
 918      */
 919     bKeepLR = (fr->bTwinRange && !bNS);
 920     for (i = 0; (i < egNR); i++)
 921     {
 922         if (!(bKeepLR && bMaster && (i == egCOULLR || i == egLJLR)))
 923         {
 924             for (j = 0; (j < enerd->grpp.nener); j++)
 925             {
 926                 enerd->grpp.ener[i][j] = 0.0;
 927             }
 928         }
 929     }
 930     for (i = 0; i < efptNR; i++)
 931     {
 932         enerd->dvdl_lin[i]    = 0.0;
 933         enerd->dvdl_nonlin[i] = 0.0;
 934     }
 935
 936     /* Normal potential energy components */
 937     for (i = 0; (i <= F_EPOT); i++)
 938     {
 939         enerd->term[i] = 0.0;
 940     }
 941     /* Initialize the dVdlambda term with the long range contribution */
 942     /* Initialize the dvdl term with the long range contribution */
 943     enerd->term[F_DVDL]            = 0.0;
 944     enerd->term[F_DVDL_COUL]       = 0.0;
 945     enerd->term[F_DVDL_VDW]        = 0.0;
 946     enerd->term[F_DVDL_BONDED]     = 0.0;
 947     enerd->term[F_DVDL_RESTRAINT]  = 0.0;
 948     enerd->term[F_DKDL]            = 0.0;
 949     if (enerd->n_lambda > 0)
 950     {
 951         for (i = 0; i < enerd->n_lambda; i++)
 952         {
 953             enerd->enerpart_lambda[i] = 0.0;
 954         }
 955     }
 956     /* reset foreign energy data - separate function since we also call it elsewhere */
 957     reset_foreign_enerdata(enerd);
 958 }