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         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 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                        gmx_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     gmx_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     set_pbc(&pbc,fr->ePBC,box);
 161
 162     /* Reset box */
 163     for(i=0; (i<DIM); i++)
 164     {
 165         box_size[i]=box[i][i];
 166     }
 167
 168     bSepDVDL=(fr->bSepDVDL && do_per_step(step,ir->nstlog));
 169     debug_gmx();
 170
 171     /* do QMMM first if requested */
 172     if(fr->bQMMM)
 173     {
 174         enerd->term[F_EQM] = calculate_QMMM(cr,x,f,fr,md);
 175     }
 176
 177     if (bSepDVDL)
 178     {
 179         fprintf(fplog,"Step %s: non-bonded V and dVdl for node %d:\n",
 180                 gmx_step_str(step,buf),cr->nodeid);
 181     }
 182
 183     /* Call the short range functions all in one go. */
 184     GMX_MPE_LOG(ev_do_fnbf_start);
 185
 186     dvdlambda = 0;
 187
 188 #ifdef GMX_MPI
 189     /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
 190 #define TAKETIME FALSE
 191     if (TAKETIME)
 192     {
 193         MPI_Barrier(cr->mpi_comm_mygroup);
 194         t0=MPI_Wtime();
 195     }
 196 #endif
 197
 198     if (ir->nwall)
 199     {
 200         dvdlambda = do_walls(ir,fr,box,md,x,f,lambda,
 201                              enerd->grpp.ener[egLJSR],nrnb);
 202         PRINT_SEPDVDL("Walls",0.0,dvdlambda);
 203         enerd->dvdl_lin += dvdlambda;
 204     }
 205
 206         /* If doing GB, reset dvda and calculate the Born radii */
 207         if (ir->implicit_solvent)
 208         {
 209                 /* wallcycle_start(wcycle,ewcGB); */
 210
 211                 for(i=0;i<born->nr;i++)
 212                 {
 213                         fr->dvda[i]=0;
 214                 }
 215
 216                 if(bBornRadii)
 217                 {
 218                         calc_gb_rad(cr,fr,ir,top,atype,x,&(fr->gblist),born,md,nrnb);
 219                 }
 220
 221                 /* wallcycle_stop(wcycle, ewcGB); */
 222         }
 223
 224     where();
 225     donb_flags = 0;
 226     if (flags & GMX_FORCE_FORCES)
 227     {
 228         donb_flags |= GMX_DONB_FORCES;
 229     }
 230     do_nonbonded(cr,fr,x,f,md,excl,
 231                  fr->bBHAM ?
 232                  enerd->grpp.ener[egBHAMSR] :
 233                  enerd->grpp.ener[egLJSR],
 234                  enerd->grpp.ener[egCOULSR],
 235                                  enerd->grpp.ener[egGB],box_size,nrnb,
 236                  lambda,&dvdlambda,-1,-1,donb_flags);
 237     /* If we do foreign lambda and we have soft-core interactions
 238      * we have to recalculate the (non-linear) energies contributions.
 239      */
 240     if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) && ir->sc_alpha != 0)
 241     {
 242         init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
 243
 244         for(i=0; i<enerd->n_lambda; i++)
 245         {
 246             lam_i = (i==0 ? lambda : ir->flambda[i-1]);
 247             dvdl_dum = 0;
 248             reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 249             do_nonbonded(cr,fr,x,f,md,excl,
 250                          fr->bBHAM ?
 251                          ed_lam.grpp.ener[egBHAMSR] :
 252                          ed_lam.grpp.ener[egLJSR],
 253                          ed_lam.grpp.ener[egCOULSR],
 254                          enerd->grpp.ener[egGB], box_size,nrnb,
 255                          lam_i,&dvdl_dum,-1,-1,
 256                          GMX_DONB_FOREIGNLAMBDA);
 257             sum_epot(&ir->opts,&ed_lam);
 258             enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 259         }
 260         destroy_enerdata(&ed_lam);
 261     }
 262     where();
 263
 264         /* If we are doing GB, calculate bonded forces and apply corrections
 265          * to the solvation forces */
 266         if (ir->implicit_solvent)  {
 267                 calc_gb_forces(cr,md,born,top,atype,x,f,fr,idef,
 268                        ir->gb_algorithm,nrnb,bBornRadii,&pbc,graph,enerd);
 269     }
 270
 271 #ifdef GMX_MPI
 272     if (TAKETIME)
 273     {
 274         t1=MPI_Wtime();
 275         fr->t_fnbf += t1-t0;
 276     }
 277 #endif
 278
 279     if (ir->sc_alpha != 0)
 280     {
 281         enerd->dvdl_nonlin += dvdlambda;
 282     }
 283     else
 284     {
 285         enerd->dvdl_lin    += dvdlambda;
 286     }
 287     Vsr = 0;
 288     if (bSepDVDL)
 289     {
 290         for(i=0; i<enerd->grpp.nener; i++)
 291         {
 292             Vsr +=
 293                 (fr->bBHAM ?
 294                  enerd->grpp.ener[egBHAMSR][i] :
 295                  enerd->grpp.ener[egLJSR][i])
 296                 + enerd->grpp.ener[egCOULSR][i] + enerd->grpp.ener[egGB][i];
 297         }
 298     }
 299     PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",Vsr,dvdlambda);
 300     debug_gmx();
 301
 302     GMX_MPE_LOG(ev_do_fnbf_finish);
 303
 304     if (debug)
 305     {
 306         pr_rvecs(debug,0,"fshift after SR",fr->fshift,SHIFTS);
 307     }
 308
 309     /* Shift the coordinates. Must be done before bonded forces and PPPM,
 310      * but is also necessary for SHAKE and update, therefore it can NOT
 311      * go when no bonded forces have to be evaluated.
 312      */
 313
 314     /* Here sometimes we would not need to shift with NBFonly,
 315      * but we do so anyhow for consistency of the returned coordinates.
 316      */
 317     if (graph)
 318     {
 319         shift_self(graph,box,x);
 320         if (TRICLINIC(box))
 321         {
 322             inc_nrnb(nrnb,eNR_SHIFTX,2*graph->nnodes);
 323         }
 324         else
 325         {
 326             inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);
 327         }
 328     }
 329     /* Check whether we need to do bondeds or correct for exclusions */
 330     if (fr->bMolPBC &&
 331         ((flags & GMX_FORCE_BONDED)
 332          || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))
 333     {
 334         /* Since all atoms are in the rectangular or triclinic unit-cell,
 335          * only single box vector shifts (2 in x) are required.
 336          */
 337         set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
 338     }
 339     debug_gmx();
 340
 341     if (flags & GMX_FORCE_BONDED)
 342     {
 343         GMX_MPE_LOG(ev_calc_bonds_start);
 344         calc_bonds(fplog,cr->ms,
 345                    idef,x,hist,f,fr,&pbc,graph,enerd,nrnb,lambda,md,fcd,
 346                    DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,
 347                    fr->bSepDVDL && do_per_step(step,ir->nstlog),step);
 348
 349         /* Check if we have to determine energy differences
 350          * at foreign lambda's.
 351          */
 352         if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) &&
 353             idef->ilsort != ilsortNO_FE)
 354         {
 355             if (idef->ilsort != ilsortFE_SORTED)
 356             {
 357                 gmx_incons("The bonded interactions are not sorted for free energy");
 358             }
 359             init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
 360
 361             for(i=0; i<enerd->n_lambda; i++)
 362             {
 363                 lam_i = (i==0 ? lambda : ir->flambda[i-1]);
 364                 dvdl_dum = 0;
 365                 reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
 366                 calc_bonds_lambda(fplog,
 367                                   idef,x,fr,&pbc,graph,&ed_lam,nrnb,lam_i,md,
 368                                   fcd,
 369                                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
 370                 sum_epot(&ir->opts,&ed_lam);
 371                 enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
 372             }
 373             destroy_enerdata(&ed_lam);
 374         }
 375         debug_gmx();
 376         GMX_MPE_LOG(ev_calc_bonds_finish);
 377     }
 378
 379     where();
 380
 381     *cycles_pme = 0;
 382     if (EEL_FULL(fr->eeltype))
 383     {
 384         bSB = (ir->nwall == 2);
 385         if (bSB)
 386         {
 387             copy_mat(box,boxs);
 388             svmul(ir->wall_ewald_zfac,boxs[ZZ],boxs[ZZ]);
 389             box_size[ZZ] *= ir->wall_ewald_zfac;
 390         }
 391
 392         clear_mat(fr->vir_el_recip);
 393
 394         if (fr->bEwald)
 395         {
 396             if (fr->n_tpi == 0)
 397             {
 398                 dvdlambda = 0;
 399                 Vcorr = ewald_LRcorrection(fplog,md->start,md->start+md->homenr,
 400                                            cr,fr,
 401                                            md->chargeA,
 402                                            md->nChargePerturbed ? md->chargeB : NULL,
 403                                            excl,x,bSB ? boxs : box,mu_tot,
 404                                            ir->ewald_geometry,
 405                                            ir->epsilon_surface,
 406                                            lambda,&dvdlambda,&vdip,&vcharge);
 407                 PRINT_SEPDVDL("Ewald excl./charge/dip. corr.",Vcorr,dvdlambda);
 408                 enerd->dvdl_lin += dvdlambda;
 409             }
 410             else
 411             {
 412                 if (ir->ewald_geometry != eewg3D || ir->epsilon_surface != 0)
 413                 {
 414                     gmx_fatal(FARGS,"TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
 415                 }
 416                 /* The TPI molecule does not have exclusions with the rest
 417                  * of the system and no intra-molecular PME grid contributions
 418                  * will be calculated in gmx_pme_calc_energy.
 419                  */
 420                 Vcorr = 0;
 421             }
 422         }
 423         else
 424         {
 425             Vcorr = shift_LRcorrection(fplog,md->start,md->homenr,cr,fr,
 426                                        md->chargeA,excl,x,TRUE,box,
 427                                        fr->vir_el_recip);
 428         }
 429
 430         dvdlambda = 0;
 431         status = 0;
 432         switch (fr->eeltype)
 433         {
 434         case eelPPPM:
 435             status = gmx_pppm_do(fplog,fr->pmedata,FALSE,x,fr->f_novirsum,
 436                                  md->chargeA,
 437                                  box_size,fr->phi,cr,md->start,md->homenr,
 438                                  nrnb,ir->pme_order,&Vlr);
 439             break;
 440         case eelPME:
 441         case eelPMESWITCH:
 442         case eelPMEUSER:
 443         case eelPMEUSERSWITCH:
 444             if (cr->duty & DUTY_PME)
 445             {
 446                 if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
 447                 {
 448                     pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
 449                     if (flags & GMX_FORCE_FORCES)
 450                     {
 451                         pme_flags |= GMX_PME_CALC_F;
 452                     }
 453                     if (flags & GMX_FORCE_VIRIAL)
 454                     {
 455                         pme_flags |= GMX_PME_CALC_ENER_VIR;
 456                     }
 457                     wallcycle_start(wcycle,ewcPMEMESH);
 458                     status = gmx_pme_do(fr->pmedata,
 459                                         md->start,md->homenr - fr->n_tpi,
 460                                         x,fr->f_novirsum,
 461                                         md->chargeA,md->chargeB,
 462                                         bSB ? boxs : box,cr,
 463                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
 464                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
 465                                         nrnb,wcycle,
 466                                         fr->vir_el_recip,fr->ewaldcoeff,
 467                                         &Vlr,lambda,&dvdlambda,
 468                                         pme_flags);
 469                     *cycles_pme = wallcycle_stop(wcycle,ewcPMEMESH);
 470
 471                     /* We should try to do as little computation after
 472                      * this as possible, because parallel PME synchronizes
 473                      * the nodes, so we want all load imbalance of the rest
 474                      * of the force calculation to be before the PME call.
 475                      * DD load balancing is done on the whole time of
 476                      * the force call (without PME).
 477                      */
 478                 }
 479                 if (fr->n_tpi > 0)
 480                 {
 481                     /* Determine the PME grid energy of the test molecule
 482                      * with the PME grid potential of the other charges.
 483                      */
 484                     gmx_pme_calc_energy(fr->pmedata,fr->n_tpi,
 485                                         x + md->homenr - fr->n_tpi,
 486                                         md->chargeA + md->homenr - fr->n_tpi,
 487                                         &Vlr);
 488                 }
 489                 PRINT_SEPDVDL("PME mesh",Vlr,dvdlambda);
 490             }
 491             else
 492             {
 493                 /* Energies and virial are obtained later from the PME nodes */
 494                 /* but values have to be zeroed out here */
 495                 Vlr=0.0;
 496             }
 497             break;
 498         case eelEWALD:
 499             Vlr = do_ewald(fplog,FALSE,ir,x,fr->f_novirsum,
 500                            md->chargeA,md->chargeB,
 501                            box_size,cr,md->homenr,
 502                            fr->vir_el_recip,fr->ewaldcoeff,
 503                            lambda,&dvdlambda,fr->ewald_table);
 504             PRINT_SEPDVDL("Ewald long-range",Vlr,dvdlambda);
 505             break;
 506         default:
 507             Vlr = 0;
 508             gmx_fatal(FARGS,"No such electrostatics method implemented %s",
 509                       eel_names[fr->eeltype]);
 510         }
 511         if (status != 0)
 512         {
 513             gmx_fatal(FARGS,"Error %d in long range electrostatics routine %s",
 514                       status,EELTYPE(fr->eeltype));
 515                 }
 516         enerd->dvdl_lin += dvdlambda;
 517         enerd->term[F_COUL_RECIP] = Vlr + Vcorr;
 518         if (debug)
 519         {
 520             fprintf(debug,"Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
 521                     Vlr,Vcorr,enerd->term[F_COUL_RECIP]);
 522             pr_rvecs(debug,0,"vir_el_recip after corr",fr->vir_el_recip,DIM);
 523             pr_rvecs(debug,0,"fshift after LR Corrections",fr->fshift,SHIFTS);
 524         }
 525     }
 526     else
 527     {
 528         if (EEL_RF(fr->eeltype))
 529         {
 530             dvdlambda = 0;
 531
 532             if (fr->eeltype != eelRF_NEC)
 533             {
 534                 enerd->term[F_RF_EXCL] =
 535                     RF_excl_correction(fplog,fr,graph,md,excl,x,f,
 536                                        fr->fshift,&pbc,lambda,&dvdlambda);
 537             }
 538
 539             enerd->dvdl_lin += dvdlambda;
 540             PRINT_SEPDVDL("RF exclusion correction",
 541                           enerd->term[F_RF_EXCL],dvdlambda);
 542         }
 543     }
 544     where();
 545     debug_gmx();
 546
 547     if (debug)
 548     {
 549         print_nrnb(debug,nrnb);
 550     }
 551     debug_gmx();
 552
 553 #ifdef GMX_MPI
 554     if (TAKETIME)
 555     {
 556         t2=MPI_Wtime();
 557         MPI_Barrier(cr->mpi_comm_mygroup);
 558         t3=MPI_Wtime();
 559         fr->t_wait += t3-t2;
 560         if (fr->timesteps == 11)
 561         {
 562             fprintf(stderr,"* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
 563                     cr->nodeid, gmx_step_str(fr->timesteps,buf),
 564                     100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
 565                     (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);
 566         }
 567         fr->timesteps++;
 568     }
 569 #endif
 570
 571     if (debug)
 572     {
 573         pr_rvecs(debug,0,"fshift after bondeds",fr->fshift,SHIFTS);
 574     }
 575
 576     GMX_MPE_LOG(ev_force_finish);
 577
 578 }
 579
 580 void init_enerdata(int ngener,int n_flambda,gmx_enerdata_t *enerd)
 581 {
 582     int i,n2;
 583
 584     for(i=0; i<F_NRE; i++)
 585     {
 586         enerd->term[i] = 0;
 587     }
 588
 589     n2=ngener*ngener;
 590     if (debug)
 591     {
 592         fprintf(debug,"Creating %d sized group matrix for energies\n",n2);
 593     }
 594     enerd->grpp.nener = n2;
 595     for(i=0; (i<egNR); i++)
 596     {
 597         snew(enerd->grpp.ener[i],n2);
 598     }
 599
 600     if (n_flambda)
 601     {
 602         enerd->n_lambda = 1 + n_flambda;
 603         snew(enerd->enerpart_lambda,enerd->n_lambda);
 604     }
 605     else
 606     {
 607         enerd->n_lambda = 0;
 608     }
 609 }
 610
 611 void destroy_enerdata(gmx_enerdata_t *enerd)
 612 {
 613     int i;
 614
 615     for(i=0; (i<egNR); i++)
 616     {
 617         sfree(enerd->grpp.ener[i]);
 618     }
 619
 620     if (enerd->n_lambda)
 621     {
 622         sfree(enerd->enerpart_lambda);
 623     }
 624 }
 625
 626 static real sum_v(int n,real v[])
 627 {
 628   real t;
 629   int  i;
 630
 631   t = 0.0;
 632   for(i=0; (i<n); i++)
 633     t = t + v[i];
 634
 635   return t;
 636 }
 637
 638 void sum_epot(t_grpopts *opts,gmx_enerdata_t *enerd)
 639 {
 640   gmx_grppairener_t *grpp;
 641   real *epot;
 642   int i;
 643
 644   grpp = &enerd->grpp;
 645   epot = enerd->term;
 646
 647   /* Accumulate energies */
 648   epot[F_COUL_SR]  = sum_v(grpp->nener,grpp->ener[egCOULSR]);
 649   epot[F_LJ]       = sum_v(grpp->nener,grpp->ener[egLJSR]);
 650   epot[F_LJ14]     = sum_v(grpp->nener,grpp->ener[egLJ14]);
 651   epot[F_COUL14]   = sum_v(grpp->nener,grpp->ener[egCOUL14]);
 652   epot[F_COUL_LR]  = sum_v(grpp->nener,grpp->ener[egCOULLR]);
 653   epot[F_LJ_LR]    = sum_v(grpp->nener,grpp->ener[egLJLR]);
 654   /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
 655   epot[F_GBPOL]   += sum_v(grpp->nener,grpp->ener[egGB]);
 656
 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,gmx_bool bNS,
 709                     gmx_enerdata_t *enerd,
 710                     gmx_bool bMaster)
 711 {
 712   gmx_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 }