src/mdlib/force.c

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