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