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 "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "macros.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "physics.h"
  48 #include "force.h"
  49 #include "nonbonded.h"
  50 #include "names.h"
  51 #include "network.h"
  52 #include "pbc.h"
  53 #include "ns.h"
  54 #include "nrnb.h"
  55 #include "bondf.h"
  56 #include "mshift.h"
  57 #include "txtdump.h"
  58 #include "coulomb.h"
  59 #include "pppm.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
  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         bool       bFillGrid,
  82         bool       bDoLongRange,
  83         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 void do_force_lowlevel(FILE       *fplog,   gmx_large_int_t step,
 117                        t_forcerec *fr,      t_inputrec *ir,
 118                        t_idef     *idef,    t_commrec  *cr,
 119                        t_nrnb     *nrnb,    gmx_wallcycle_t wcycle,
 120                        t_mdatoms  *md,
 121                        t_grpopts  *opts,
 122                        rvec       x[],      history_t  *hist,
 123                        rvec       f[],
 124                        gmx_enerdata_t *enerd,
 125                        t_fcdata   *fcd,
 126                        gmx_mtop_t     *mtop,
 127                        gmx_localtop_t *top,
 128                        gmx_genborn_t *born,
 129                        t_atomtypes *atype,
 130                        bool       bBornRadii,
 131                        matrix     box,
 132                        real       lambda,
 133                        t_graph    *graph,
 134                        t_blocka   *excl,
 135                        rvec       mu_tot[],
 136                        int        flags,
 137                        float      *cycles_pme)
 138 {
 139     int     i,status;
 140     int     donb_flags;
 141     bool    bDoEpot,bSepDVDL,bSB;
 142     int     pme_flags;
 143     matrix  boxs;
 144     rvec    box_size;
 145     real    dvdlambda,Vsr,Vlr,Vcorr=0,vdip,vcharge;
 146     t_pbc   pbc;
 147     real    dvdgb;
 148     char    buf[22];
 149     gmx_enerdata_t ed_lam;
 150     double  lam_i;
 151     real    dvdl_dum;
 152
 153 #ifdef GMX_MPI
 154     double  t0=0.0,t1,t2,t3; /* time measurement for coarse load balancing */
 155 #endif
 156
 157 #define PRINT_SEPDVDL(s,v,dvdl) if (bSepDVDL) fprintf(fplog,sepdvdlformat,s,v,dvdl);
 158
 159     GMX_MPE_LOG(ev_force_start);
 160
 161     /* Reset box */
 162     for(i=0; (i<DIM); i++)
 163     {
 164         box_size[i]=box[i][i];
 165     }
 166
 167     bSepDVDL=(fr->bSepDVDL && do_per_step(step,ir->nstlog));
 168     debug_gmx();
 169
 170     /* do QMMM first if requested */
 171     if(fr->bQMMM)
 172     {
 173         enerd->term[F_EQM] = calculate_QMMM(cr,x,f,fr,md);
 174     }
 175
 176     if (bSepDVDL)
 177     {
 178         fprintf(fplog,"Step %s: non-bonded V and dVdl for node %d:\n",
 179                 gmx_step_str(step,buf),cr->nodeid);
 180     }
 181
 182     /* Call the short range functions all in one go. */
 183     GMX_MPE_LOG(ev_do_fnbf_start);
 184
 185     dvdlambda = 0;
 186
 187 #ifdef GMX_MPI
 188     /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
 189 #define TAKETIME FALSE
 190     if (TAKETIME)
 191     {
 192         MPI_Barrier(cr->mpi_comm_mygroup);
 193         t0=MPI_Wtime();
 194     }
 195 #endif
 196
 197     if (ir->nwall)
 198     {
 199         dvdlambda = do_walls(ir,fr,box,md,x,f,lambda,
 200                              enerd->grpp.ener[egLJSR],nrnb);
 201         PRINT_SEPDVDL("Walls",0.0,dvdlambda);
 202         enerd->dvdl_lin += dvdlambda;
 203     }
 204
 205         /* If doing GB, reset dvda and calculate the Born radii */
 206         if (ir->implicit_solvent)
 207         {
 208                 /* wallcycle_start(wcycle,ewcGB); */
 209
 210                 for(i=0;i<born->nr;i++)
 211                 {
 212                         fr->dvda[i]=0;
 213                 }
 214
 215                 if(bBornRadii)
 216                 {
 217                         calc_gb_rad(cr,fr,ir,top,atype,x,&(fr->gblist),born,md,nrnb);
 218                 }
 219
 220                 /* wallcycle_stop(wcycle, ewcGB); */
 221         }
 222
 223     where();
 224     donb_flags = 0;
 225     if (flags & GMX_FORCE_FORCES)
 226     {
 227         donb_flags |= GMX_DONB_FORCES;
 228     }
 229     do_nonbonded(cr,fr,x,f,md,excl,
 230                  fr->bBHAM ?
 231                  enerd->grpp.ener[egBHAMSR] :
 232                  enerd->grpp.ener[egLJSR],
 233                  enerd->grpp.ener[egCOULSR],
 234                                  enerd->grpp.ener[egGB],box_size,nrnb,
 235                  lambda,&dvdlambda,-1,-1,donb_flags);
 236     /* If we do foreign lambda and we have soft-core interactions
 237      * we have to recalculate the (non-linear) energies contributions.
 238      */
 239     if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) && ir->sc_alpha != 0)
 240     {
 241         init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
 242
 243         for(i=0; i<enerd->n_lambda; i++)
 244         {
 245             lam_i = (i==0 ? lambda : ir->flambda[i-1]);
 246             dvdl_dum = 0;
 247             reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 248             do_nonbonded(cr,fr,x,f,md,excl,
 249                          fr->bBHAM ?
 250                          ed_lam.grpp.ener[egBHAMSR] :
 251                          ed_lam.grpp.ener[egLJSR],
 252                          ed_lam.grpp.ener[egCOULSR],
 253                          enerd->grpp.ener[egGB], box_size,nrnb,
 254                          lam_i,&dvdl_dum,-1,-1,
 255                          GMX_DONB_FOREIGNLAMBDA);
 256             sum_epot(&ir->opts,&ed_lam);
 257             enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 258         }
 259         destroy_enerdata(&ed_lam);
 260     }
 261     where();
 262
 263         /* If we are doing GB, calculate bonded forces and apply corrections
 264          * to the solvation forces */
 265         if (ir->implicit_solvent)  {
 266                 dvdgb = calc_gb_forces(cr,md,born,top,atype,x,f,fr,idef,
 267                                                            ir->gb_algorithm,nrnb,bBornRadii,&pbc,graph);
 268
 269                 enerd->term[F_GB12]+=dvdgb;
 270
 271                 /* Also add the nonbonded GB potential energy (only from one energy group currently) */
 272                 enerd->term[F_GB12]+=enerd->grpp.ener[egGB][0];
 273         }
 274
 275 #ifdef GMX_MPI
 276     if (TAKETIME)
 277     {
 278         t1=MPI_Wtime();
 279         fr->t_fnbf += t1-t0;
 280     }
 281 #endif
 282
 283     if (ir->sc_alpha != 0)
 284     {
 285         enerd->dvdl_nonlin += dvdlambda;
 286     }
 287     else
 288     {
 289         enerd->dvdl_lin    += dvdlambda;
 290     }
 291     Vsr = 0;
 292     if (bSepDVDL)
 293     {
 294         for(i=0; i<enerd->grpp.nener; i++)
 295         {
 296             Vsr +=
 297                 (fr->bBHAM ?
 298                  enerd->grpp.ener[egBHAMSR][i] :
 299                  enerd->grpp.ener[egLJSR][i])
 300                 + enerd->grpp.ener[egCOULSR][i];
 301         }
 302     }
 303     PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",Vsr,dvdlambda);
 304     debug_gmx();
 305
 306     GMX_MPE_LOG(ev_do_fnbf_finish);
 307
 308     if (debug)
 309     {
 310         pr_rvecs(debug,0,"fshift after SR",fr->fshift,SHIFTS);
 311     }
 312
 313     /* Shift the coordinates. Must be done before bonded forces and PPPM,
 314      * but is also necessary for SHAKE and update, therefore it can NOT
 315      * go when no bonded forces have to be evaluated.
 316      */
 317
 318     /* Here sometimes we would not need to shift with NBFonly,
 319      * but we do so anyhow for consistency of the returned coordinates.
 320      */
 321     if (graph)
 322     {
 323         shift_self(graph,box,x);
 324         if (TRICLINIC(box))
 325         {
 326             inc_nrnb(nrnb,eNR_SHIFTX,2*graph->nnodes);
 327         }
 328         else
 329         {
 330             inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);
 331         }
 332     }
 333     /* Check whether we need to do bondeds or correct for exclusions */
 334     if (fr->bMolPBC &&
 335         ((flags & GMX_FORCE_BONDED)
 336          || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))
 337     {
 338         /* Since all atoms are in the rectangular or triclinic unit-cell,
 339          * only single box vector shifts (2 in x) are required.
 340          */
 341         set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
 342     }
 343     debug_gmx();
 344
 345     if (flags & GMX_FORCE_BONDED)
 346     {
 347         GMX_MPE_LOG(ev_calc_bonds_start);
 348         calc_bonds(fplog,cr->ms,
 349                    idef,x,hist,f,fr,&pbc,graph,enerd,nrnb,lambda,md,fcd,
 350                    DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,
 351                    fr->bSepDVDL && do_per_step(step,ir->nstlog),step);
 352
 353         /* Check if we have to determine energy differences
 354          * at foreign lambda's.
 355          */
 356         if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) &&
 357             idef->ilsort != ilsortNO_FE)
 358         {
 359             if (idef->ilsort != ilsortFE_SORTED)
 360             {
 361                 gmx_incons("The bonded interactions are not sorted for free energy");
 362             }
 363             init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
 364
 365             for(i=0; i<enerd->n_lambda; i++)
 366             {
 367                 lam_i = (i==0 ? lambda : ir->flambda[i-1]);
 368                 dvdl_dum = 0;
 369                 reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 370                 calc_bonds_lambda(fplog,
 371                                   idef,x,fr,&pbc,graph,&ed_lam,nrnb,lam_i,md,
 372                                   fcd,
 373                                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
 374                 sum_epot(&ir->opts,&ed_lam);
 375                 enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 376             }
 377             destroy_enerdata(&ed_lam);
 378         }
 379         debug_gmx();
 380         GMX_MPE_LOG(ev_calc_bonds_finish);
 381     }
 382
 383     where();
 384
 385     *cycles_pme = 0;
 386     if (EEL_FULL(fr->eeltype))
 387     {
 388         bSB = (ir->nwall == 2);
 389         if (bSB)
 390         {
 391             copy_mat(box,boxs);
 392             svmul(ir->wall_ewald_zfac,boxs[ZZ],boxs[ZZ]);
 393             box_size[ZZ] *= ir->wall_ewald_zfac;
 394         }
 395
 396         clear_mat(fr->vir_el_recip);
 397
 398         if (fr->bEwald)
 399         {
 400             if (fr->n_tpi == 0)
 401             {
 402                 dvdlambda = 0;
 403                 Vcorr = ewald_LRcorrection(fplog,md->start,md->start+md->homenr,
 404                                            cr,fr,
 405                                            md->chargeA,
 406                                            md->nChargePerturbed ? md->chargeB : NULL,
 407                                            excl,x,bSB ? boxs : box,mu_tot,
 408                                            ir->ewald_geometry,
 409                                            ir->epsilon_surface,
 410                                            lambda,&dvdlambda,&vdip,&vcharge);
 411                 PRINT_SEPDVDL("Ewald excl./charge/dip. corr.",Vcorr,dvdlambda);
 412                 enerd->dvdl_lin += dvdlambda;
 413             }
 414             else
 415             {
 416                 if (ir->ewald_geometry != eewg3D || ir->epsilon_surface != 0)
 417                 {
 418                     gmx_fatal(FARGS,"TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
 419                 }
 420                 /* The TPI molecule does not have exclusions with the rest
 421                  * of the system and no intra-molecular PME grid contributions
 422                  * will be calculated in gmx_pme_calc_energy.
 423                  */
 424                 Vcorr = 0;
 425             }
 426         }
 427         else
 428         {
 429             Vcorr = shift_LRcorrection(fplog,md->start,md->homenr,cr,fr,
 430                                        md->chargeA,excl,x,TRUE,box,
 431                                        fr->vir_el_recip);
 432         }
 433
 434         dvdlambda = 0;
 435         status = 0;
 436         switch (fr->eeltype)
 437         {
 438         case eelPPPM:
 439             status = gmx_pppm_do(fplog,fr->pmedata,FALSE,x,fr->f_novirsum,
 440                                  md->chargeA,
 441                                  box_size,fr->phi,cr,md->start,md->homenr,
 442                                  nrnb,ir->pme_order,&Vlr);
 443             break;
 444         case eelPME:
 445         case eelPMESWITCH:
 446         case eelPMEUSER:
 447             if (cr->duty & DUTY_PME)
 448             {
 449                 if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
 450                 {
 451                     pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
 452                     if (flags & GMX_FORCE_FORCES)
 453                     {
 454                         pme_flags |= GMX_PME_CALC_F;
 455                     }
 456                     if (flags & GMX_FORCE_VIRIAL)
 457                     {
 458                         pme_flags |= GMX_PME_CALC_ENER_VIR;
 459                     }
 460                     wallcycle_start(wcycle,ewcPMEMESH);
 461                     status = gmx_pme_do(fr->pmedata,
 462                                         md->start,md->homenr - fr->n_tpi,
 463                                         x,fr->f_novirsum,
 464                                         md->chargeA,md->chargeB,
 465                                         bSB ? boxs : box,cr,
 466                                         DOMAINDECOMP(cr) ? dd_pme_maxshift0(cr->dd) : 0,
 467                                         DOMAINDECOMP(cr) ? dd_pme_maxshift1(cr->dd) : 0,
 468                                         nrnb,wcycle,
 469                                         fr->vir_el_recip,fr->ewaldcoeff,
 470                                         &Vlr,lambda,&dvdlambda,
 471                                         pme_flags);
 472                     *cycles_pme = wallcycle_stop(wcycle,ewcPMEMESH);
 473
 474                     /* We should try to do as little computation after
 475                      * this as possible, because parallel PME synchronizes
 476                      * the nodes, so we want all load imbalance of the rest
 477                      * of the force calculation to be before the PME call.
 478                      * DD load balancing is done on the whole time of
 479                      * the force call (without PME).
 480                      */
 481                 }
 482                 if (fr->n_tpi > 0)
 483                 {
 484                     /* Determine the PME grid energy of the test molecule
 485                      * with the PME grid potential of the other charges.
 486                      */
 487                     gmx_pme_calc_energy(fr->pmedata,fr->n_tpi,
 488                                         x + md->homenr - fr->n_tpi,
 489                                         md->chargeA + md->homenr - fr->n_tpi,
 490                                         &Vlr);
 491                 }
 492                 PRINT_SEPDVDL("PME mesh",Vlr,dvdlambda);
 493             }
 494             else
 495             {
 496                 /* Energies and virial are obtained later from the PME nodes */
 497                 /* but values have to be zeroed out here */
 498                 Vlr=0.0;
 499             }
 500             break;
 501         case eelEWALD:
 502             Vlr = do_ewald(fplog,FALSE,ir,x,fr->f_novirsum,
 503                            md->chargeA,md->chargeB,
 504                            box_size,cr,md->homenr,
 505                            fr->vir_el_recip,fr->ewaldcoeff,
 506                            lambda,&dvdlambda,fr->ewald_table);
 507             PRINT_SEPDVDL("Ewald long-range",Vlr,dvdlambda);
 508             break;
 509         default:
 510             Vlr = 0;
 511             gmx_fatal(FARGS,"No such electrostatics method implemented %s",
 512                       eel_names[fr->eeltype]);
 513         }
 514         if (status != 0)
 515         {
 516             gmx_fatal(FARGS,"Error %d in long range electrostatics routine %s",
 517                       status,EELTYPE(fr->eeltype));
 518                 }
 519         enerd->dvdl_lin += dvdlambda;
 520         enerd->term[F_COUL_RECIP] = Vlr + Vcorr;
 521         if (debug)
 522         {
 523             fprintf(debug,"Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
 524                     Vlr,Vcorr,enerd->term[F_COUL_RECIP]);
 525             pr_rvecs(debug,0,"vir_el_recip after corr",fr->vir_el_recip,DIM);
 526             pr_rvecs(debug,0,"fshift after LR Corrections",fr->fshift,SHIFTS);
 527         }
 528     }
 529     else
 530     {
 531         if (EEL_RF(fr->eeltype))
 532         {
 533             dvdlambda = 0;
 534
 535             if (fr->eeltype != eelRF_NEC)
 536             {
 537                 enerd->term[F_RF_EXCL] =
 538                     RF_excl_correction(fplog,fr,graph,md,excl,x,f,
 539                                        fr->fshift,&pbc,lambda,&dvdlambda);
 540             }
 541
 542             enerd->dvdl_lin += dvdlambda;
 543             PRINT_SEPDVDL("RF exclusion correction",
 544                           enerd->term[F_RF_EXCL],dvdlambda);
 545         }
 546     }
 547     where();
 548     debug_gmx();
 549
 550     if (debug)
 551     {
 552         print_nrnb(debug,nrnb);
 553     }
 554     debug_gmx();
 555
 556 #ifdef GMX_MPI
 557     if (TAKETIME)
 558     {
 559         t2=MPI_Wtime();
 560         MPI_Barrier(cr->mpi_comm_mygroup);
 561         t3=MPI_Wtime();
 562         fr->t_wait += t3-t2;
 563         if (fr->timesteps == 11)
 564         {
 565             fprintf(stderr,"* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
 566                     cr->nodeid, gmx_step_str(fr->timesteps,buf),
 567                     100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
 568                     (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);
 569         }
 570         fr->timesteps++;
 571     }
 572 #endif
 573
 574     if (debug)
 575     {
 576         pr_rvecs(debug,0,"fshift after bondeds",fr->fshift,SHIFTS);
 577     }
 578
 579     GMX_MPE_LOG(ev_force_finish);
 580
 581 }
 582
 583 void init_enerdata(int ngener,int n_flambda,gmx_enerdata_t *enerd)
 584 {
 585     int i,n2;
 586
 587     for(i=0; i<F_NRE; i++)
 588     {
 589         enerd->term[i] = 0;
 590     }
 591
 592     n2=ngener*ngener;
 593     if (debug)
 594     {
 595         fprintf(debug,"Creating %d sized group matrix for energies\n",n2);
 596     }
 597     enerd->grpp.nener = n2;
 598     for(i=0; (i<egNR); i++)
 599     {
 600         snew(enerd->grpp.ener[i],n2);
 601     }
 602
 603     if (n_flambda)
 604     {
 605         enerd->n_lambda = 1 + n_flambda;
 606         snew(enerd->enerpart_lambda,enerd->n_lambda);
 607     }
 608     else
 609     {
 610         enerd->n_lambda = 0;
 611     }
 612 }
 613
 614 void destroy_enerdata(gmx_enerdata_t *enerd)
 615 {
 616     int i;
 617
 618     for(i=0; (i<egNR); i++)
 619     {
 620         sfree(enerd->grpp.ener[i]);
 621     }
 622
 623     if (enerd->n_lambda)
 624     {
 625         sfree(enerd->enerpart_lambda);
 626     }
 627 }
 628
 629 static real sum_v(int n,real v[])
 630 {
 631   real t;
 632   int  i;
 633
 634   t = 0.0;
 635   for(i=0; (i<n); i++)
 636     t = t + v[i];
 637
 638   return t;
 639 }
 640
 641 void sum_epot(t_grpopts *opts,gmx_enerdata_t *enerd)
 642 {
 643   gmx_grppairener_t *grpp;
 644   real *epot;
 645   int i;
 646
 647   grpp = &enerd->grpp;
 648   epot = enerd->term;
 649
 650   /* Accumulate energies */
 651   epot[F_COUL_SR]  = sum_v(grpp->nener,grpp->ener[egCOULSR]);
 652   epot[F_LJ]       = sum_v(grpp->nener,grpp->ener[egLJSR]);
 653   epot[F_LJ14]     = sum_v(grpp->nener,grpp->ener[egLJ14]);
 654   epot[F_COUL14]   = sum_v(grpp->nener,grpp->ener[egCOUL14]);
 655   epot[F_COUL_LR]  = sum_v(grpp->nener,grpp->ener[egCOULLR]);
 656   epot[F_LJ_LR]    = sum_v(grpp->nener,grpp->ener[egLJLR]);
 657 /* lattice part of LR doesnt belong to any group
 658  * and has been added earlier
 659  */
 660   epot[F_BHAM]     = sum_v(grpp->nener,grpp->ener[egBHAMSR]);
 661   epot[F_BHAM_LR]  = sum_v(grpp->nener,grpp->ener[egBHAMLR]);
 662
 663   epot[F_EPOT] = 0;
 664   for(i=0; (i<F_EPOT); i++)
 665     if (i != F_DISRESVIOL && i != F_ORIRESDEV && i != F_DIHRESVIOL)
 666       epot[F_EPOT] += epot[i];
 667 }
 668
 669 void sum_dhdl(gmx_enerdata_t *enerd,double lambda,t_inputrec *ir)
 670 {
 671     int i;
 672     double dlam,dhdl_lin;
 673
 674     enerd->term[F_DVDL] = enerd->dvdl_lin + enerd->dvdl_nonlin;
 675
 676     if (debug)
 677     {
 678         fprintf(debug,"dvdl: %f, non-linear %f + linear %f\n",
 679                 enerd->term[F_DVDL],enerd->dvdl_nonlin,enerd->dvdl_lin);
 680     }
 681
 682     /* Notes on the foreign lambda free energy difference evaluation:
 683      * Adding the potential and ekin terms that depend linearly on lambda
 684      * as delta lam * dvdl to the energy differences is exact.
 685      * For the constraint dvdl this is not exact, but we have no other option.
 686      * For the non-bonded LR term we assume that the soft-core (if present)
 687      * no longer affects the energy beyond the short-range cut-off,
 688      * which is a very good approximation (except for exotic settings).
 689      */
 690     for(i=1; i<enerd->n_lambda; i++)
 691     {
 692         dlam = (ir->flambda[i-1] - lambda);
 693         dhdl_lin =
 694             enerd->dvdl_lin + enerd->term[F_DKDL] + enerd->term[F_DHDL_CON];
 695         if (debug)
 696         {
 697             fprintf(debug,"enerdiff lam %g: non-linear %f linear %f*%f\n",
 698                     ir->flambda[i-1],
 699                     enerd->enerpart_lambda[i] - enerd->enerpart_lambda[0],
 700                     dlam,dhdl_lin);
 701         }
 702         enerd->enerpart_lambda[i] += dlam*dhdl_lin;
 703
 704     }
 705 }
 706
 707 void reset_enerdata(t_grpopts *opts,
 708                     t_forcerec *fr,bool bNS,
 709                     gmx_enerdata_t *enerd,
 710                     bool bMaster)
 711 {
 712   bool bKeepLR;
 713   int  i,j;
 714
 715   /* First reset all energy components, except for the long range terms
 716    * on the master at non neighbor search steps, since the long range
 717    * terms have already been summed at the last neighbor search step.
 718    */
 719   bKeepLR = (fr->bTwinRange && !bNS);
 720   for(i=0; (i<egNR); i++) {
 721     if (!(bKeepLR && bMaster && (i == egCOULLR || i == egLJLR))) {
 722       for(j=0; (j<enerd->grpp.nener); j++)
 723         enerd->grpp.ener[i][j] = 0.0;
 724     }
 725   }
 726   enerd->dvdl_lin    = 0.0;
 727   enerd->dvdl_nonlin = 0.0;
 728
 729   /* Normal potential energy components */
 730   for(i=0; (i<=F_EPOT); i++) {
 731     enerd->term[i] = 0.0;
 732   }
 733   /* Initialize the dVdlambda term with the long range contribution */
 734   enerd->term[F_DVDL]     = 0.0;
 735   enerd->term[F_DKDL]     = 0.0;
 736   enerd->term[F_DHDL_CON] = 0.0;
 737   if (enerd->n_lambda > 0)
 738   {
 739       for(i=0; i<enerd->n_lambda; i++)
 740       {
 741           enerd->enerpart_lambda[i] = 0.0;
 742       }
 743   }
 744 }