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