src/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 "mpelogging.h"
  66 #include "gmx_omp_nthreads.h"
  67
  68 void ns(FILE *fp,
  69         t_forcerec *fr,
  70         rvec       x[],
  71         matrix     box,
  72         gmx_groups_t *groups,
  73         t_grpopts  *opts,
  74         gmx_localtop_t *top,
  75         t_mdatoms  *md,
  76         t_commrec  *cr,
  77         t_nrnb     *nrnb,
  78         real       *lambda,
  79         real       *dvdlambda,
  80         gmx_grppairener_t *grppener,
  81         gmx_bool       bFillGrid,
  82         gmx_bool       bDoLongRange,
  83         gmx_bool       bDoForces,
  84         rvec       *f)
  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,bDoLongRange,
 101                                 bDoForces,f);
 102   if (debug)
 103     fprintf(debug,"nsearch = %d\n",nsearch);
 104
 105   /* Check whether we have to do dynamic load balancing */
 106   /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
 107     count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
 108     &(top->idef),opts->ngener);
 109   */
 110   if (fr->ns.dump_nl > 0)
 111     dump_nblist(fp,cr,fr,fr->ns.dump_nl);
 112
 113   GMX_MPE_LOG(ev_ns_finish);
 114 }
 115
 116 static void reduce_thread_forces(int n,rvec *f,
 117                                  tensor vir,
 118                                  real *Vcorr,
 119                                  int efpt_ind,real *dvdl,
 120                                  int nthreads,f_thread_t *f_t)
 121 {
 122     int t,i;
 123
 124     /* This reduction can run over any number of threads */
 125 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
 126     for(i=0; i<n; i++)
 127     {
 128         for(t=1; t<nthreads; t++)
 129         {
 130             rvec_inc(f[i],f_t[t].f[i]);
 131         }
 132     }
 133     for(t=1; t<nthreads; t++)
 134     {
 135         *Vcorr += f_t[t].Vcorr;
 136         *dvdl  += f_t[t].dvdl[efpt_ind];
 137         m_add(vir,f_t[t].vir,vir);
 138     }
 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                        t_grpopts  *opts,
 147                        rvec       x[],      history_t  *hist,
 148                        rvec       f[],
 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     gmx_enerdata_t ed_lam;
 176     double  clam_i,vlam_i;
 177     real    dvdl_dum[efptNR], dvdl, dvdl_nb[efptNR], lam_i[efptNR];
 178     real    dvdlsum;
 179
 180 #ifdef GMX_MPI
 181     double  t0=0.0,t1,t2,t3; /* time measurement for coarse load balancing */
 182 #endif
 183
 184 #define PRINT_SEPDVDL(s,v,dvdlambda) if (bSepDVDL) fprintf(fplog,sepdvdlformat,s,v,dvdlambda);
 185
 186     GMX_MPE_LOG(ev_force_start);
 187     set_pbc(&pbc,fr->ePBC,box);
 188
 189     /* reset free energy components */
 190     for (i=0;i<efptNR;i++)
 191     {
 192         dvdl_nb[i]  = 0;
 193         dvdl_dum[i] = 0;
 194     }
 195
 196     /* Reset box */
 197     for(i=0; (i<DIM); i++)
 198     {
 199         box_size[i]=box[i][i];
 200     }
 201
 202     bSepDVDL=(fr->bSepDVDL && do_per_step(step,ir->nstlog));
 203     debug_gmx();
 204
 205     /* do QMMM first if requested */
 206     if(fr->bQMMM)
 207     {
 208         enerd->term[F_EQM] = calculate_QMMM(cr,x,f,fr,md);
 209     }
 210
 211     if (bSepDVDL)
 212     {
 213         fprintf(fplog,"Step %s: non-bonded V and dVdl for node %d:\n",
 214                 gmx_step_str(step,buf),cr->nodeid);
 215     }
 216
 217     /* Call the short range functions all in one go. */
 218     GMX_MPE_LOG(ev_do_fnbf_start);
 219
 220 #ifdef GMX_MPI
 221     /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
 222 #define TAKETIME FALSE
 223     if (TAKETIME)
 224     {
 225         MPI_Barrier(cr->mpi_comm_mygroup);
 226         t0=MPI_Wtime();
 227     }
 228 #endif
 229
 230     if (ir->nwall)
 231     {
 232         /* foreign lambda component for walls */
 233         dvdl = do_walls(ir,fr,box,md,x,f,lambda[efptVDW],
 234                         enerd->grpp.ener[egLJSR],nrnb);
 235         PRINT_SEPDVDL("Walls",0.0,dvdl);
 236         enerd->dvdl_lin[efptVDW] += dvdl;
 237     }
 238
 239         /* If doing GB, reset dvda and calculate the Born radii */
 240         if (ir->implicit_solvent)
 241         {
 242         wallcycle_sub_start(wcycle, ewcsNONBONDED);
 243
 244                 for(i=0;i<born->nr;i++)
 245                 {
 246                         fr->dvda[i]=0;
 247                 }
 248
 249                 if(bBornRadii)
 250                 {
 251                         calc_gb_rad(cr,fr,ir,top,atype,x,&(fr->gblist),born,md,nrnb);
 252                 }
 253
 254         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
 255         }
 256
 257     where();
 258     if (flags & GMX_FORCE_NONBONDED)
 259     {
 260         donb_flags = 0;
 261         if (flags & GMX_FORCE_FORCES)
 262         {
 263             donb_flags |= GMX_DONB_FORCES;
 264         }
 265
 266         wallcycle_sub_start(wcycle, ewcsNONBONDED);
 267         do_nonbonded(cr,fr,x,f,md,excl,
 268                     fr->bBHAM ?
 269                     enerd->grpp.ener[egBHAMSR] :
 270                     enerd->grpp.ener[egLJSR],
 271                     enerd->grpp.ener[egCOULSR],
 272                     enerd->grpp.ener[egGB],box_size,nrnb,
 273                     lambda,dvdl_nb,-1,-1,donb_flags);
 274         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
 275     }
 276
 277     /* If we do foreign lambda and we have soft-core interactions
 278      * we have to recalculate the (non-linear) energies contributions.
 279      */
 280     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL) && fepvals->sc_alpha != 0)
 281     {
 282         wallcycle_sub_start(wcycle, ewcsNONBONDED);
 283         init_enerdata(mtop->groups.grps[egcENER].nr,fepvals->n_lambda,&ed_lam);
 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_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 292             do_nonbonded(cr,fr,x,f,md,excl,
 293                          fr->bBHAM ?
 294                          ed_lam.grpp.ener[egBHAMSR] :
 295                          ed_lam.grpp.ener[egLJSR],
 296                          ed_lam.grpp.ener[egCOULSR],
 297                          enerd->grpp.ener[egGB], box_size,nrnb,
 298                          lam_i,dvdl_dum,-1,-1,
 299                          GMX_DONB_FOREIGNLAMBDA);
 300             sum_epot(&ir->opts,&ed_lam);
 301             enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 302         }
 303         destroy_enerdata(&ed_lam);
 304         wallcycle_sub_stop(wcycle, ewcsNONBONDED);
 305     }
 306     where();
 307
 308         /* If we are doing GB, calculate bonded forces and apply corrections
 309          * to the solvation forces */
 310     /* MRS: Eventually, many need to include free energy contribution here! */
 311         if (ir->implicit_solvent)
 312     {
 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             init_enerdata(mtop->groups.grps[egcENER].nr,fepvals->n_lambda,&ed_lam);
 424
 425             for(i=0; i<enerd->n_lambda; i++)
 426             {
 427                 reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 428                 for (j=0;j<efptNR;j++)
 429                 {
 430                     lam_i[j] = (i==0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 431                 }
 432                 calc_bonds_lambda(fplog,idef,x,fr,&pbc,graph,&ed_lam,nrnb,lam_i,md,
 433                                   fcd,DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
 434                 sum_epot(&ir->opts,&ed_lam);
 435                 enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 436             }
 437             destroy_enerdata(&ed_lam);
 438         }
 439         debug_gmx();
 440         GMX_MPE_LOG(ev_calc_bonds_finish);
 441         wallcycle_sub_stop(wcycle, ewcsBONDED);
 442     }
 443
 444     where();
 445
 446     *cycles_pme = 0;
 447     if (EEL_FULL(fr->eeltype))
 448     {
 449         bSB = (ir->nwall == 2);
 450         if (bSB)
 451         {
 452             copy_mat(box,boxs);
 453             svmul(ir->wall_ewald_zfac,boxs[ZZ],boxs[ZZ]);
 454             box_size[ZZ] *= ir->wall_ewald_zfac;
 455         }
 456
 457         clear_mat(fr->vir_el_recip);
 458
 459         if (fr->bEwald)
 460         {
 461             Vcorr = 0;
 462             dvdl  = 0;
 463
 464             /* With the Verlet scheme exclusion forces are calculated
 465              * in the non-bonded kernel.
 466              */
 467             /* The TPI molecule does not have exclusions with the rest
 468              * of the system and no intra-molecular PME grid contributions
 469              * will be calculated in gmx_pme_calc_energy.
 470              */
 471             if ((ir->cutoff_scheme == ecutsGROUP && fr->n_tpi == 0) ||
 472                 ir->ewald_geometry != eewg3D ||
 473                 ir->epsilon_surface != 0)
 474             {
 475                 int nthreads,t;
 476
 477                 wallcycle_sub_start(wcycle, ewcsEWALD_CORRECTION);
 478
 479                 if (fr->n_tpi > 0)
 480                 {
 481                     gmx_fatal(FARGS,"TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
 482                 }
 483
 484                 nthreads = gmx_omp_nthreads_get(emntBonded);
 485 #pragma omp parallel for num_threads(nthreads) schedule(static)
 486                 for(t=0; t<nthreads; t++)
 487                 {
 488                     int s,e,i;
 489                     rvec *fnv;
 490                     tensor *vir;
 491                     real *Vcorrt,*dvdlt;
 492                     if (t == 0)
 493                     {
 494                         fnv    = fr->f_novirsum;
 495                         vir    = &fr->vir_el_recip;
 496                         Vcorrt = &Vcorr;
 497                         dvdlt  = &dvdl;
 498                     }
 499                     else
 500                     {
 501                         fnv    = fr->f_t[t].f;
 502                         vir    = &fr->f_t[t].vir;
 503                         Vcorrt = &fr->f_t[t].Vcorr;
 504                         dvdlt  = &fr->f_t[t].dvdl[efptCOUL];
 505                         for(i=0; i<fr->natoms_force; i++)
 506                         {
 507                             clear_rvec(fnv[i]);
 508                         }
 509                         clear_mat(*vir);
 510                     }
 511                     *dvdlt = 0;
 512                     *Vcorrt =
 513                         ewald_LRcorrection(fplog,
 514                                            fr->excl_load[t],fr->excl_load[t+1],
 515                                            cr,t,fr,
 516                                            md->chargeA,
 517                                            md->nChargePerturbed ? md->chargeB : NULL,
 518                                            ir->cutoff_scheme != ecutsVERLET,
 519                                            excl,x,bSB ? boxs : box,mu_tot,
 520                                            ir->ewald_geometry,
 521                                            ir->epsilon_surface,
 522                                            fnv,*vir,
 523                                            lambda[efptCOUL],dvdlt);
 524                 }
 525                 if (nthreads > 1)
 526                 {
 527                     reduce_thread_forces(fr->natoms_force,fr->f_novirsum,
 528                                          fr->vir_el_recip,
 529                                          &Vcorr,efptCOUL,&dvdl,
 530                                          nthreads,fr->f_t);
 531                 }
 532
 533                 wallcycle_sub_stop(wcycle, ewcsEWALD_CORRECTION);
 534             }
 535
 536             if (fr->n_tpi == 0)
 537             {
 538                 Vcorr += ewald_charge_correction(cr,fr,lambda[efptCOUL],box,
 539                                                  &dvdl,fr->vir_el_recip);
 540             }
 541
 542             PRINT_SEPDVDL("Ewald excl./charge/dip. corr.",Vcorr,dvdl);
 543             enerd->dvdl_lin[efptCOUL] += dvdl;
 544         }
 545
 546         status = 0;
 547         dvdl = 0;
 548         switch (fr->eeltype)
 549         {
 550         case eelPME:
 551         case eelPMESWITCH:
 552         case eelPMEUSER:
 553         case eelPMEUSERSWITCH:
 554         case eelP3M_AD:
 555             if (cr->duty & DUTY_PME)
 556             {
 557                 assert(fr->n_tpi >= 0);
 558                 if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
 559                 {
 560                     pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
 561                     if (flags & GMX_FORCE_FORCES)
 562                     {
 563                         pme_flags |= GMX_PME_CALC_F;
 564                     }
 565                     if (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY))
 566                     {
 567                         pme_flags |= GMX_PME_CALC_ENER_VIR;
 568                     }
 569                     if (fr->n_tpi > 0)
 570                     {
 571                         /* We don't calculate f, but we do want the potential */
 572                         pme_flags |= GMX_PME_CALC_POT;
 573                     }
 574                     wallcycle_start(wcycle,ewcPMEMESH);
 575                     status = gmx_pme_do(fr->pmedata,
 576                                         md->start,md->homenr - fr->n_tpi,
 577                                         x,fr->f_novirsum,
 578                                         md->chargeA,md->chargeB,
 579                                         bSB ? boxs : box,cr,
 580                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
 581                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
 582                                         nrnb,wcycle,
 583                                         fr->vir_el_recip,fr->ewaldcoeff,
 584                                         &Vlr,lambda[efptCOUL],&dvdl,
 585                                         pme_flags);
 586                     *cycles_pme = wallcycle_stop(wcycle,ewcPMEMESH);
 587
 588                     /* We should try to do as little computation after
 589                      * this as possible, because parallel PME synchronizes
 590                      * the nodes, so we want all load imbalance of the rest
 591                      * of the force calculation to be before the PME call.
 592                      * DD load balancing is done on the whole time of
 593                      * the force call (without PME).
 594                      */
 595                 }
 596                 if (fr->n_tpi > 0)
 597                 {
 598                     /* Determine the PME grid energy of the test molecule
 599                      * with the PME grid potential of the other charges.
 600                      */
 601                     gmx_pme_calc_energy(fr->pmedata,fr->n_tpi,
 602                                         x + md->homenr - fr->n_tpi,
 603                                         md->chargeA + md->homenr - fr->n_tpi,
 604                                         &Vlr);
 605                 }
 606                 PRINT_SEPDVDL("PME mesh",Vlr,dvdl);
 607             }
 608             else
 609             {
 610                 /* Energies and virial are obtained later from the PME nodes */
 611                 /* but values have to be zeroed out here */
 612                 Vlr=0.0;
 613             }
 614             break;
 615         case eelEWALD:
 616             Vlr = do_ewald(fplog,FALSE,ir,x,fr->f_novirsum,
 617                            md->chargeA,md->chargeB,
 618                            box_size,cr,md->homenr,
 619                            fr->vir_el_recip,fr->ewaldcoeff,
 620                            lambda[efptCOUL],&dvdl,fr->ewald_table);
 621             PRINT_SEPDVDL("Ewald long-range",Vlr,dvdl);
 622             break;
 623         default:
 624             Vlr = 0;
 625             gmx_fatal(FARGS,"No such electrostatics method implemented %s",
 626                       eel_names[fr->eeltype]);
 627         }
 628         if (status != 0)
 629         {
 630             gmx_fatal(FARGS,"Error %d in long range electrostatics routine %s",
 631                       status,EELTYPE(fr->eeltype));
 632                 }
 633         enerd->dvdl_lin[efptCOUL] += dvdl;
 634         enerd->term[F_COUL_RECIP] = Vlr + Vcorr;
 635         if (debug)
 636         {
 637             fprintf(debug,"Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
 638                     Vlr,Vcorr,enerd->term[F_COUL_RECIP]);
 639             pr_rvecs(debug,0,"vir_el_recip after corr",fr->vir_el_recip,DIM);
 640             pr_rvecs(debug,0,"fshift after LR Corrections",fr->fshift,SHIFTS);
 641         }
 642     }
 643     else
 644     {
 645         if (EEL_RF(fr->eeltype))
 646         {
 647             /* With the Verlet scheme exclusion forces are calculated
 648              * in the non-bonded kernel.
 649              */
 650             if (ir->cutoff_scheme != ecutsVERLET && fr->eeltype != eelRF_NEC)
 651             {
 652                 dvdl = 0;
 653                 enerd->term[F_RF_EXCL] =
 654                     RF_excl_correction(fplog,fr,graph,md,excl,x,f,
 655                                        fr->fshift,&pbc,lambda[efptCOUL],&dvdl);
 656             }
 657
 658             enerd->dvdl_lin[efptCOUL] += dvdl;
 659             PRINT_SEPDVDL("RF exclusion correction",
 660                           enerd->term[F_RF_EXCL],dvdl);
 661         }
 662     }
 663     where();
 664     debug_gmx();
 665
 666     if (debug)
 667     {
 668         print_nrnb(debug,nrnb);
 669     }
 670     debug_gmx();
 671
 672 #ifdef GMX_MPI
 673     if (TAKETIME)
 674     {
 675         t2=MPI_Wtime();
 676         MPI_Barrier(cr->mpi_comm_mygroup);
 677         t3=MPI_Wtime();
 678         fr->t_wait += t3-t2;
 679         if (fr->timesteps == 11)
 680         {
 681             fprintf(stderr,"* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
 682                     cr->nodeid, gmx_step_str(fr->timesteps,buf),
 683                     100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
 684                     (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);
 685         }
 686         fr->timesteps++;
 687     }
 688 #endif
 689
 690     if (debug)
 691     {
 692         pr_rvecs(debug,0,"fshift after bondeds",fr->fshift,SHIFTS);
 693     }
 694
 695     GMX_MPE_LOG(ev_force_finish);
 696
 697 }
 698
 699 void init_enerdata(int ngener,int n_lambda,gmx_enerdata_t *enerd)
 700 {
 701     int i,n2;
 702
 703     for(i=0; i<F_NRE; i++)
 704     {
 705         enerd->term[i] = 0;
 706     }
 707
 708
 709     for(i=0; i<efptNR; i++) {
 710         enerd->dvdl_lin[i]  = 0;
 711         enerd->dvdl_nonlin[i]  = 0;
 712     }
 713
 714     n2=ngener*ngener;
 715     if (debug)
 716     {
 717         fprintf(debug,"Creating %d sized group matrix for energies\n",n2);
 718     }
 719     enerd->grpp.nener = n2;
 720     for(i=0; (i<egNR); i++)
 721     {
 722         snew(enerd->grpp.ener[i],n2);
 723     }
 724
 725     if (n_lambda)
 726     {
 727         enerd->n_lambda = 1 + n_lambda;
 728         snew(enerd->enerpart_lambda,enerd->n_lambda);
 729     }
 730     else
 731     {
 732         enerd->n_lambda = 0;
 733     }
 734 }
 735
 736 void destroy_enerdata(gmx_enerdata_t *enerd)
 737 {
 738     int i;
 739
 740     for(i=0; (i<egNR); i++)
 741     {
 742         sfree(enerd->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_enerdata_t *enerd)
 764 {
 765   gmx_grppairener_t *grpp;
 766   real *epot;
 767   int i;
 768
 769   grpp = &enerd->grpp;
 770   epot = enerd->term;
 771
 772   /* Accumulate energies */
 773   epot[F_COUL_SR]  = sum_v(grpp->nener,grpp->ener[egCOULSR]);
 774   epot[F_LJ]       = sum_v(grpp->nener,grpp->ener[egLJSR]);
 775   epot[F_LJ14]     = sum_v(grpp->nener,grpp->ener[egLJ14]);
 776   epot[F_COUL14]   = sum_v(grpp->nener,grpp->ener[egCOUL14]);
 777   epot[F_COUL_LR]  = sum_v(grpp->nener,grpp->ener[egCOULLR]);
 778   epot[F_LJ_LR]    = sum_v(grpp->nener,grpp->ener[egLJLR]);
 779   /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
 780   epot[F_GBPOL]   += sum_v(grpp->nener,grpp->ener[egGB]);
 781
 782 /* lattice part of LR doesnt belong to any group
 783  * and has been added earlier
 784  */
 785   epot[F_BHAM]     = sum_v(grpp->nener,grpp->ener[egBHAMSR]);
 786   epot[F_BHAM_LR]  = sum_v(grpp->nener,grpp->ener[egBHAMLR]);
 787
 788   epot[F_EPOT] = 0;
 789   for(i=0; (i<F_EPOT); i++)
 790   {
 791       if (i != F_DISRESVIOL && i != F_ORIRESDEV)
 792       {
 793           epot[F_EPOT] += epot[i];
 794       }
 795   }
 796 }
 797
 798 void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
 799 {
 800     int i,j,index;
 801     double dlam;
 802
 803     enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW];  /* include dispersion correction */
 804     enerd->term[F_DVDL] = 0.0;
 805     for (i=0;i<efptNR;i++)
 806     {
 807         if (fepvals->separate_dvdl[i])
 808         {
 809             /* could this be done more readably/compactly? */
 810             switch (i) {
 811             case (efptCOUL):
 812                 index = F_DVDL_COUL;
 813                 break;
 814             case (efptVDW):
 815                 index = F_DVDL_VDW;
 816                 break;
 817             case (efptBONDED):
 818                 index = F_DVDL_BONDED;
 819                 break;
 820             case (efptRESTRAINT):
 821                 index = F_DVDL_RESTRAINT;
 822                 break;
 823             case (efptMASS):
 824                 index = F_DKDL;
 825                 break;
 826             default:
 827                 index = F_DVDL;
 828                 break;
 829             }
 830             enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 831             if (debug)
 832             {
 833                 fprintf(debug,"dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
 834                         efpt_names[i],i,enerd->term[index],enerd->dvdl_nonlin[i],enerd->dvdl_lin[i]);
 835             }
 836         }
 837         else
 838         {
 839             enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 840             if (debug)
 841             {
 842                 fprintf(debug,"dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
 843                         efpt_names[0],i,enerd->term[F_DVDL],enerd->dvdl_nonlin[i],enerd->dvdl_lin[i]);
 844             }
 845         }
 846     }
 847
 848     /* Notes on the foreign lambda free energy difference evaluation:
 849      * Adding the potential and ekin terms that depend linearly on lambda
 850      * as delta lam * dvdl to the energy differences is exact.
 851      * For the constraints this is not exact, but we have no other option
 852      * without literally changing the lengths and reevaluating the energies at each step.
 853      * (try to remedy this post 4.6 - MRS)
 854      * For the non-bonded LR term we assume that the soft-core (if present)
 855      * no longer affects the energy beyond the short-range cut-off,
 856      * which is a very good approximation (except for exotic settings).
 857      * (investigate how to overcome this post 4.6 - MRS)
 858      */
 859
 860     for(i=0; i<fepvals->n_lambda; i++)
 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 contributions to the current lambda, because
 867            it's automatically zero */
 868
 869         /* first kinetic energy term */
 870         dlam = (fepvals->all_lambda[efptMASS][i] - lambda[efptMASS]);
 871
 872         enerd->enerpart_lambda[i+1] += enerd->term[F_DKDL]*dlam;
 873
 874         for (j=0;j<efptNR;j++)
 875         {
 876             if (j==efptMASS) {continue;} /* no other mass term to worry about */
 877
 878             dlam = (fepvals->all_lambda[j][i]-lambda[j]);
 879             enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
 880             if (debug)
 881             {
 882                 fprintf(debug,"enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
 883                         fepvals->all_lambda[j][i],efpt_names[j],
 884                         (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),
 885                         dlam,enerd->dvdl_lin[j]);
 886             }
 887         }
 888     }
 889 }
 890
 891 void reset_enerdata(t_grpopts *opts,
 892                     t_forcerec *fr,gmx_bool bNS,
 893                     gmx_enerdata_t *enerd,
 894                     gmx_bool bMaster)
 895 {
 896     gmx_bool bKeepLR;
 897     int  i,j;
 898
 899     /* First reset all energy components, except for the long range terms
 900      * on the master at non neighbor search steps, since the long range
 901      * terms have already been summed at the last neighbor search step.
 902      */
 903     bKeepLR = (fr->bTwinRange && !bNS);
 904     for(i=0; (i<egNR); i++) {
 905         if (!(bKeepLR && bMaster && (i == egCOULLR || i == egLJLR))) {
 906             for(j=0; (j<enerd->grpp.nener); j++)
 907                 enerd->grpp.ener[i][j] = 0.0;
 908         }
 909     }
 910     for (i=0;i<efptNR;i++)
 911     {
 912         enerd->dvdl_lin[i]    = 0.0;
 913         enerd->dvdl_nonlin[i] = 0.0;
 914     }
 915
 916     /* Normal potential energy components */
 917     for(i=0; (i<=F_EPOT); i++) {
 918         enerd->term[i] = 0.0;
 919     }
 920     /* Initialize the dVdlambda term with the long range contribution */
 921     /* Initialize the dvdl term with the long range contribution */
 922     enerd->term[F_DVDL]            = 0.0;
 923     enerd->term[F_DVDL_COUL]       = 0.0;
 924     enerd->term[F_DVDL_VDW]        = 0.0;
 925     enerd->term[F_DVDL_BONDED]     = 0.0;
 926     enerd->term[F_DVDL_RESTRAINT]  = 0.0;
 927     enerd->term[F_DKDL]            = 0.0;
 928     if (enerd->n_lambda > 0)
 929     {
 930         for(i=0; i<enerd->n_lambda; i++)
 931         {
 932             enerd->enerpart_lambda[i] = 0.0;
 933         }
 934     }
 935 }