src/gromacs/gmxlib/nonbonded/nonbonded.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  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #ifdef GMX_THREAD_MPI
  41 #include <thread_mpi.h>
  42 #endif
  43
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include "typedefs.h"
  47 #include "txtdump.h"
  48 #include "smalloc.h"
  49 #include "ns.h"
  50 #include "vec.h"
  51 #include "maths.h"
  52 #include "macros.h"
  53 #include "string2.h"
  54 #include "force.h"
  55 #include "names.h"
  56 #include "main.h"
  57 #include "xvgr.h"
  58 #include "gmx_fatal.h"
  59 #include "physics.h"
  60 #include "force.h"
  61 #include "bondf.h"
  62 #include "nrnb.h"
  63 #include "smalloc.h"
  64 #include "nonbonded.h"
  65
  66 #include "nb_kernel.h"
  67 #include "nb_free_energy.h"
  68 #include "nb_generic.h"
  69 #include "nb_generic_cg.h"
  70 #include "nb_generic_adress.h"
  71
  72 /* Different default (c) and accelerated interaction-specific kernels */
  73 #include "nb_kernel_c/nb_kernel_c.h"
  74
  75 #if (defined GMX_CPU_ACCELERATION_X86_SSE2) && !(defined GMX_DOUBLE)
  76 #    include "nb_kernel_sse2_single/nb_kernel_sse2_single.h"
  77 #endif
  78 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1) && !(defined GMX_DOUBLE)
  79 #    include "nb_kernel_sse4_1_single/nb_kernel_sse4_1_single.h"
  80 #endif
  81 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA) && !(defined GMX_DOUBLE)
  82 #    include "nb_kernel_avx_128_fma_single/nb_kernel_avx_128_fma_single.h"
  83 #endif
  84 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256) && !(defined GMX_DOUBLE)
  85 #    include "nb_kernel_avx_256_single/nb_kernel_avx_256_single.h"
  86 #endif
  87 #if (defined GMX_CPU_ACCELERATION_X86_SSE2 && defined GMX_DOUBLE)
  88 #    include "nb_kernel_sse2_double/nb_kernel_sse2_double.h"
  89 #endif
  90 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1 && defined GMX_DOUBLE)
  91 #    include "nb_kernel_sse4_1_double/nb_kernel_sse4_1_double.h"
  92 #endif
  93 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA && defined GMX_DOUBLE)
  94 #    include "nb_kernel_avx_128_fma_double/nb_kernel_avx_128_fma_double.h"
  95 #endif
  96 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256 && defined GMX_DOUBLE)
  97 #    include "nb_kernel_avx_256_double/nb_kernel_avx_256_double.h"
  98 #endif
  99
 100
 101 #ifdef GMX_THREAD_MPI
 102 static tMPI_Thread_mutex_t nonbonded_setup_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
 103 #endif
 104 static gmx_bool            nonbonded_setup_done  = FALSE;
 105
 106
 107 void
 108 gmx_nonbonded_setup(FILE *         fplog,
 109                     t_forcerec *   fr,
 110                     gmx_bool       bGenericKernelOnly)
 111 {
 112 #ifdef GMX_THREAD_MPI
 113     tMPI_Thread_mutex_lock(&nonbonded_setup_mutex);
 114 #endif
 115     /* Here we are guaranteed only one thread made it. */
 116     if(nonbonded_setup_done==FALSE)
 117     {
 118         if(bGenericKernelOnly==FALSE)
 119         {
 120             /* Add the generic kernels to the structure stored statically in nb_kernel.c */
 121             nb_kernel_list_add_kernels(kernellist_c,kernellist_c_size);
 122
 123             if(!(fr!=NULL && fr->use_cpu_acceleration==FALSE))
 124             {
 125                 /* Add interaction-specific kernels for different architectures */
 126                 /* Single precision */
 127 #if (defined GMX_CPU_ACCELERATION_X86_SSE2) && !(defined GMX_DOUBLE)
 128                 nb_kernel_list_add_kernels(kernellist_sse2_single,kernellist_sse2_single_size);
 129 #endif
 130 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1) && !(defined GMX_DOUBLE)
 131                 nb_kernel_list_add_kernels(kernellist_sse4_1_single,kernellist_sse4_1_single_size);
 132 #endif
 133 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA) && !(defined GMX_DOUBLE)
 134                 nb_kernel_list_add_kernels(kernellist_avx_128_fma_single,kernellist_avx_128_fma_single_size);
 135 #endif
 136 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256) && !(defined GMX_DOUBLE)
 137                 nb_kernel_list_add_kernels(kernellist_avx_256_single,kernellist_avx_256_single_size);
 138 #endif
 139                 /* Double precision */
 140 #if (defined GMX_CPU_ACCELERATION_X86_SSE2 && defined GMX_DOUBLE)
 141                 nb_kernel_list_add_kernels(kernellist_sse2_double,kernellist_sse2_double_size);
 142 #endif
 143 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1 && defined GMX_DOUBLE)
 144                 nb_kernel_list_add_kernels(kernellist_sse4_1_double,kernellist_sse4_1_double_size);
 145 #endif
 146 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA && defined GMX_DOUBLE)
 147                 nb_kernel_list_add_kernels(kernellist_avx_128_fma_double,kernellist_avx_128_fma_double_size);
 148 #endif
 149 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256 && defined GMX_DOUBLE)
 150                 nb_kernel_list_add_kernels(kernellist_avx_256_double,kernellist_avx_256_double_size);
 151 #endif
 152                 ; /* empty statement to avoid a completely empty block */
 153             }
 154         }
 155         /* Create a hash for faster lookups */
 156         nb_kernel_list_hash_init();
 157
 158         nonbonded_setup_done=TRUE;
 159     }
 160 #ifdef GMX_THREAD_MPI
 161     tMPI_Thread_mutex_unlock(&nonbonded_setup_mutex);
 162 #endif
 163 }
 164
 165
 166
 167 void
 168 gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl)
 169 {
 170     const char *     elec;
 171     const char *     elec_mod;
 172     const char *     vdw;
 173     const char *     vdw_mod;
 174     const char *     geom;
 175     const char *     other;
 176     const char *     vf;
 177
 178     struct
 179     {
 180         const char *  arch;
 181         int           simd_padding_width;
 182     }
 183     arch_and_padding[] =
 184     {
 185         /* Single precision */
 186 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256) && !(defined GMX_DOUBLE)
 187         { "avx_256_single", 8 },
 188 #endif
 189 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA) && !(defined GMX_DOUBLE)
 190         { "avx_128_fma_single", 4 },
 191 #endif
 192 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1) && !(defined GMX_DOUBLE)
 193         { "sse4_1_single", 4 },
 194 #endif
 195 #if (defined GMX_CPU_ACCELERATION_X86_SSE2) && !(defined GMX_DOUBLE)
 196         { "sse2_single", 4 },
 197 #endif
 198         /* Double precision */
 199 #if (defined GMX_CPU_ACCELERATION_X86_AVX_256 && defined GMX_DOUBLE)
 200         { "avx_256_double", 4 },
 201 #endif
 202 #if (defined GMX_CPU_ACCELERATION_X86_AVX_128_FMA && defined GMX_DOUBLE)
 203         /* Sic. Double precision 2-way SIMD does not require neighbor list padding,
 204          * since the kernels execute a loop unrolled a factor 2, followed by
 205          * a possible single odd-element epilogue.
 206          */
 207         { "avx_128_fma_double", 1 },
 208 #endif
 209 #if (defined GMX_CPU_ACCELERATION_X86_SSE2 && defined GMX_DOUBLE)
 210         /* No padding - see comment above */
 211         { "sse2_double", 1 },
 212 #endif
 213 #if (defined GMX_CPU_ACCELERATION_X86_SSE4_1 && defined GMX_DOUBLE)
 214         /* No padding - see comment above */
 215         { "sse4_1_double", 1 },
 216 #endif
 217         { "c", 1 },
 218     };
 219     int              narch = asize(arch_and_padding);
 220     int              i;
 221
 222     if(nonbonded_setup_done==FALSE)
 223     {
 224         /* We typically call this setup routine before starting timers,
 225          * but if that has not been done for whatever reason we do it now.
 226          */
 227         gmx_nonbonded_setup(NULL,NULL,FALSE);
 228     }
 229
 230     /* Not used yet */
 231     other="";
 232
 233     nl->kernelptr_vf = NULL;
 234     nl->kernelptr_v  = NULL;
 235     nl->kernelptr_f  = NULL;
 236
 237     elec     = gmx_nbkernel_elec_names[nl->ielec];
 238     elec_mod = eintmod_names[nl->ielecmod];
 239     vdw      = gmx_nbkernel_vdw_names[nl->ivdw];
 240     vdw_mod  = eintmod_names[nl->ivdwmod];
 241     geom     = gmx_nblist_geometry_names[nl->igeometry];
 242
 243     if(nl->free_energy)
 244     {
 245         nl->kernelptr_vf = gmx_nb_free_energy_kernel;
 246         nl->kernelptr_f  = gmx_nb_free_energy_kernel;
 247         nl->simd_padding_width = 1;
 248     }
 249     else if(!gmx_strcasecmp_min(geom,"CG-CG"))
 250     {
 251         nl->kernelptr_vf = gmx_nb_generic_cg_kernel;
 252         nl->kernelptr_f  = gmx_nb_generic_cg_kernel;
 253         nl->simd_padding_width = 1;
 254     }
 255     else
 256     {
 257         /* Try to find a specific kernel first */
 258
 259         for(i=0;i<narch && nl->kernelptr_vf==NULL ;i++)
 260         {
 261             nl->kernelptr_vf = nb_kernel_list_findkernel(log,arch_and_padding[i].arch,elec,elec_mod,vdw,vdw_mod,geom,other,"PotentialAndForce");
 262             nl->simd_padding_width = arch_and_padding[i].simd_padding_width;
 263         }
 264         for(i=0;i<narch && nl->kernelptr_f==NULL ;i++)
 265         {
 266             nl->kernelptr_f = nb_kernel_list_findkernel(log,arch_and_padding[i].arch,elec,elec_mod,vdw,vdw_mod,geom,other,"Force");
 267             nl->simd_padding_width = arch_and_padding[i].simd_padding_width;
 268
 269             /* If there is not force-only optimized kernel, is there a potential & force one? */
 270             if(nl->kernelptr_f == NULL)
 271             {
 272                 nl->kernelptr_f  = nb_kernel_list_findkernel(NULL,arch_and_padding[i].arch,elec,elec_mod,vdw,vdw_mod,geom,other,"PotentialAndForce");
 273                 nl->simd_padding_width = arch_and_padding[i].simd_padding_width;
 274             }
 275         }
 276
 277         /* Give up, pick a generic one instead */
 278         if(nl->kernelptr_vf==NULL)
 279         {
 280             nl->kernelptr_vf = gmx_nb_generic_kernel;
 281             nl->kernelptr_f  = gmx_nb_generic_kernel;
 282             nl->simd_padding_width = 1;
 283             if(debug)
 284             {
 285                 fprintf(debug,
 286                         "WARNING - Slow generic NB kernel used for neighborlist with\n"
 287                         "    Elec: '%s', Modifier: '%s'\n"
 288                         "    Vdw:  '%s', Modifier: '%s'\n"
 289                         "    Geom: '%s', Other: '%s'\n\n",
 290                         elec,elec_mod,vdw,vdw_mod,geom,other);
 291             }
 292         }
 293     }
 294
 295     return;
 296 }
 297
 298 void do_nonbonded(t_commrec *cr,t_forcerec *fr,
 299                   rvec x[],rvec f_shortrange[],rvec f_longrange[],t_mdatoms *mdatoms,t_blocka *excl,
 300                   gmx_grppairener_t *grppener,rvec box_size,
 301                   t_nrnb *nrnb,real *lambda, real *dvdl,
 302                   int nls,int eNL,int flags)
 303 {
 304         t_nblist *        nlist;
 305         int               n,n0,n1,i,i0,i1,sz,range;
 306         t_nblists *       nblists;
 307     nb_kernel_data_t  kernel_data;
 308     nb_kernel_t *     kernelptr=NULL;
 309     rvec *            f;
 310
 311     kernel_data.flags                   = flags;
 312     kernel_data.exclusions              = excl;
 313     kernel_data.lambda                  = lambda;
 314     kernel_data.dvdl                    = dvdl;
 315
 316     if(fr->bAllvsAll)
 317     {
 318         return;
 319     }
 320
 321     if (eNL >= 0)
 322     {
 323                 i0 = eNL;
 324                 i1 = i0+1;
 325     }
 326     else
 327     {
 328                 i0 = 0;
 329                 i1 = eNL_NR;
 330         }
 331
 332         if (nls >= 0)
 333         {
 334                 n0 = nls;
 335                 n1 = nls+1;
 336         }
 337         else
 338         {
 339                 n0 = 0;
 340                 n1 = fr->nnblists;
 341         }
 342
 343         for(n=n0; (n<n1); n++)
 344         {
 345                 nblists = &fr->nblists[n];
 346
 347         kernel_data.table_elec              = &nblists->table_elec;
 348         kernel_data.table_vdw               = &nblists->table_vdw;
 349         kernel_data.table_elec_vdw          = &nblists->table_elec_vdw;
 350
 351         for(range=0;range<2;range++)
 352         {
 353             /* Are we doing short/long-range? */
 354             if(range==0)
 355             {
 356                 /* Short-range */
 357                 if(!(flags & GMX_NONBONDED_DO_SR))
 358                 {
 359                     continue;
 360                 }
 361                 kernel_data.energygrp_elec          = grppener->ener[egCOULSR];
 362                 kernel_data.energygrp_vdw           = grppener->ener[fr->bBHAM ? egBHAMSR : egLJSR];
 363                 kernel_data.energygrp_polarization  = grppener->ener[egGB];
 364                 nlist = nblists->nlist_sr;
 365                 f                                   = f_shortrange;
 366             }
 367             else if(range==1)
 368             {
 369                 /* Long-range */
 370                 if(!(flags & GMX_NONBONDED_DO_LR))
 371                 {
 372                     continue;
 373                 }
 374                 kernel_data.energygrp_elec          = grppener->ener[egCOULLR];
 375                 kernel_data.energygrp_vdw           = grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR];
 376                 kernel_data.energygrp_polarization  = grppener->ener[egGB];
 377                 nlist = nblists->nlist_lr;
 378                 f                                   = f_longrange;
 379             }
 380
 381             for(i=i0; (i<i1); i++)
 382             {
 383                 if (nlist[i].nri > 0)
 384                 {
 385                     if(flags & GMX_NONBONDED_DO_POTENTIAL)
 386                     {
 387                         /* Potential and force */
 388                         kernelptr = (nb_kernel_t *)nlist[i].kernelptr_vf;
 389                     }
 390                     else
 391                     {
 392                         /* Force only, no potential */
 393                         kernelptr = (nb_kernel_t *)nlist[i].kernelptr_f;
 394                     }
 395
 396                     if(nlist[i].free_energy==0 && (flags & GMX_NONBONDED_DO_FOREIGNLAMBDA))
 397                     {
 398                         /* We don't need the non-perturbed interactions */
 399                         continue;
 400                     }
 401                     (*kernelptr)(&(nlist[i]),x,f,fr,mdatoms,&kernel_data,nrnb);
 402                  }
 403             }
 404         }
 405     }
 406 }
 407
 408 static void
 409 nb_listed_warning_rlimit(const rvec *x,int ai, int aj,int * global_atom_index,real r, real rlimit)
 410 {
 411     gmx_warning("Listed nonbonded interaction between particles %d and %d\n"
 412                 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
 413                 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
 414                 "a smaller molecule you are decoupling during a free energy calculation.\n"
 415                 "Since interactions at distances beyond the table cannot be computed,\n"
 416                 "they are skipped until they are inside the table limit again. You will\n"
 417                 "only see this message once, even if it occurs for several interactions.\n\n"
 418                 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
 419                 "probably something wrong with your system. Only change the table-extension\n"
 420                 "distance in the mdp file if you are really sure that is the reason.\n",
 421                 glatnr(global_atom_index,ai),glatnr(global_atom_index,aj),r,rlimit);
 422
 423     if (debug)
 424     {
 425         fprintf(debug,
 426                 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
 427                 x[ai][XX],x[ai][YY],x[ai][ZZ],x[aj][XX],x[aj][YY],x[aj][ZZ],
 428                 glatnr(global_atom_index,ai),glatnr(global_atom_index,aj),r);
 429     }
 430 }
 431
 432
 433
 434 /* This might logically belong better in the nb_generic.c module, but it is only
 435  * used in do_nonbonded_listed(), and we want it to be inlined there to avoid an
 436  * extra functional call for every single pair listed in the topology.
 437  */
 438 static real
 439 nb_evaluate_single(real r2, real tabscale,real *vftab,
 440                    real qq, real c6, real c12, real *velec, real *vvdw)
 441 {
 442     real       rinv,r,rtab,eps,eps2,Y,F,Geps,Heps2,Fp,VVe,FFe,VVd,FFd,VVr,FFr,fscal;
 443     int        ntab;
 444
 445     /* Do the tabulated interactions - first table lookup */
 446     rinv             = gmx_invsqrt(r2);
 447     r                = r2*rinv;
 448     rtab             = r*tabscale;
 449     ntab             = rtab;
 450     eps              = rtab-ntab;
 451     eps2             = eps*eps;
 452     ntab             = 12*ntab;
 453     /* Electrostatics */
 454     Y                = vftab[ntab];
 455     F                = vftab[ntab+1];
 456     Geps             = eps*vftab[ntab+2];
 457     Heps2            = eps2*vftab[ntab+3];
 458     Fp               = F+Geps+Heps2;
 459     VVe              = Y+eps*Fp;
 460     FFe              = Fp+Geps+2.0*Heps2;
 461     /* Dispersion */
 462     Y                = vftab[ntab+4];
 463     F                = vftab[ntab+5];
 464     Geps             = eps*vftab[ntab+6];
 465     Heps2            = eps2*vftab[ntab+7];
 466     Fp               = F+Geps+Heps2;
 467     VVd              = Y+eps*Fp;
 468     FFd              = Fp+Geps+2.0*Heps2;
 469     /* Repulsion */
 470     Y                = vftab[ntab+8];
 471     F                = vftab[ntab+9];
 472     Geps             = eps*vftab[ntab+10];
 473     Heps2            = eps2*vftab[ntab+11];
 474     Fp               = F+Geps+Heps2;
 475     VVr              = Y+eps*Fp;
 476     FFr              = Fp+Geps+2.0*Heps2;
 477
 478     *velec           = qq*VVe;
 479     *vvdw            = c6*VVd+c12*VVr;
 480
 481     fscal            = -(qq*FFe+c6*FFd+c12*FFr)*tabscale*rinv;
 482
 483     return fscal;
 484 }
 485
 486
 487 real
 488 do_nonbonded_listed(int ftype,int nbonds,
 489                 const t_iatom iatoms[],const t_iparams iparams[],
 490                 const rvec x[],rvec f[],rvec fshift[],
 491                 const t_pbc *pbc,const t_graph *g,
 492                 real *lambda, real *dvdl,
 493                 const t_mdatoms *md,
 494                 const t_forcerec *fr,gmx_grppairener_t *grppener,
 495                 int *global_atom_index)
 496 {
 497     int              ielec,ivdw;
 498     real             qq,c6,c12;
 499     rvec             dx;
 500     ivec             dt;
 501     int              i,j,itype,ai,aj,gid;
 502     int              fshift_index;
 503     real             r2,rinv;
 504     real             fscal,velec,vvdw;
 505     real *           energygrp_elec;
 506     real *           energygrp_vdw;
 507     static gmx_bool  warned_rlimit=FALSE;
 508     /* Free energy stuff */
 509     gmx_bool         bFreeEnergy;
 510     real             LFC[2],LFV[2],DLF[2],lfac_coul[2],lfac_vdw[2],dlfac_coul[2],dlfac_vdw[2];
 511     real             qqB,c6B,c12B,sigma2_def,sigma2_min;
 512
 513
 514     switch (ftype) {
 515         case F_LJ14:
 516         case F_LJC14_Q:
 517             energygrp_elec = grppener->ener[egCOUL14];
 518             energygrp_vdw  = grppener->ener[egLJ14];
 519             break;
 520         case F_LJC_PAIRS_NB:
 521             energygrp_elec = grppener->ener[egCOULSR];
 522             energygrp_vdw  = grppener->ener[egLJSR];
 523             break;
 524         default:
 525             energygrp_elec = NULL; /* Keep compiler happy */
 526             energygrp_vdw  = NULL; /* Keep compiler happy */
 527             gmx_fatal(FARGS,"Unknown function type %d in do_nonbonded14",ftype);
 528             break;
 529     }
 530
 531     if(fr->efep != efepNO)
 532     {
 533         /* Lambda factor for state A=1-lambda and B=lambda */
 534         LFC[0] = 1.0 - lambda[efptCOUL];
 535         LFV[0] = 1.0 - lambda[efptVDW];
 536         LFC[1] = lambda[efptCOUL];
 537         LFV[1] = lambda[efptVDW];
 538
 539         /*derivative of the lambda factor for state A and B */
 540         DLF[0] = -1;
 541         DLF[1] = 1;
 542
 543         /* precalculate */
 544         sigma2_def = pow(fr->sc_sigma6_def,1.0/3.0);
 545         sigma2_min = pow(fr->sc_sigma6_min,1.0/3.0);
 546
 547         for (i=0;i<2;i++)
 548         {
 549             lfac_coul[i]  = (fr->sc_power==2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 550             dlfac_coul[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power==2 ? (1-LFC[i]) : 1);
 551             lfac_vdw[i]   = (fr->sc_power==2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 552             dlfac_vdw[i]  = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power==2 ? (1-LFV[i]) : 1);
 553         }
 554     }
 555     else
 556     {
 557         sigma2_min = sigma2_def = 0;
 558     }
 559
 560     bFreeEnergy = FALSE;
 561     for(i=0; (i<nbonds); )
 562     {
 563         itype = iatoms[i++];
 564         ai    = iatoms[i++];
 565         aj    = iatoms[i++];
 566         gid   = GID(md->cENER[ai],md->cENER[aj],md->nenergrp);
 567
 568         /* Get parameters */
 569         switch (ftype) {
 570             case F_LJ14:
 571                 bFreeEnergy =
 572                 (fr->efep != efepNO &&
 573                  ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
 574                   iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
 575                   iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 576                 qq               = md->chargeA[ai]*md->chargeA[aj]*fr->epsfac*fr->fudgeQQ;
 577                 c6               = iparams[itype].lj14.c6A;
 578                 c12              = iparams[itype].lj14.c12A;
 579                 break;
 580             case F_LJC14_Q:
 581                 qq               = iparams[itype].ljc14.qi*iparams[itype].ljc14.qj*fr->epsfac*iparams[itype].ljc14.fqq;
 582                 c6               = iparams[itype].ljc14.c6;
 583                 c12              = iparams[itype].ljc14.c12;
 584                 break;
 585             case F_LJC_PAIRS_NB:
 586                 qq               = iparams[itype].ljcnb.qi*iparams[itype].ljcnb.qj*fr->epsfac;
 587                 c6               = iparams[itype].ljcnb.c6;
 588                 c12              = iparams[itype].ljcnb.c12;
 589                 break;
 590             default:
 591                 /* Cannot happen since we called gmx_fatal() above in this case */
 592                 qq = c6 = c12 = 0; /* Keep compiler happy */
 593                 break;
 594         }
 595
 596         /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
 597          * included in the general nfbp array now. This means the tables are scaled down by the
 598          * same factor, so when we use the original c6/c12 parameters from iparams[] they must
 599          * be scaled up.
 600          */
 601         c6  *= 6.0;
 602         c12 *= 12.0;
 603
 604         /* Do we need to apply full periodic boundary conditions? */
 605         if(fr->bMolPBC==TRUE)
 606         {
 607             fshift_index = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
 608         }
 609         else
 610         {
 611             fshift_index = CENTRAL;
 612             rvec_sub(x[ai],x[aj],dx);
 613         }
 614         r2           = norm2(dx);
 615
 616         if(r2>=fr->tab14.r*fr->tab14.r)
 617         {
 618             if(warned_rlimit==FALSE)
 619             {
 620                 nb_listed_warning_rlimit(x,ai,aj,global_atom_index,sqrt(r2),fr->tab14.r);
 621                 warned_rlimit=TRUE;
 622             }
 623             continue;
 624         }
 625
 626         if (bFreeEnergy)
 627         {
 628             /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
 629             qqB              = md->chargeB[ai]*md->chargeB[aj]*fr->epsfac*fr->fudgeQQ;
 630             c6B              = iparams[itype].lj14.c6B*6.0;
 631             c12B             = iparams[itype].lj14.c12B*12.0;
 632
 633             fscal            = nb_free_energy_evaluate_single(r2,fr->sc_r_power,fr->sc_alphacoul,fr->sc_alphavdw,
 634                                                               fr->tab14.scale,fr->tab14.data,qq,c6,c12,qqB,c6B,c12B,
 635                                                               LFC,LFV,DLF,lfac_coul,lfac_vdw,dlfac_coul,dlfac_vdw,
 636                                                               fr->sc_sigma6_def,fr->sc_sigma6_min,sigma2_def,sigma2_min,&velec,&vvdw,dvdl);
 637         }
 638         else
 639         {
 640             /* Evaluate tabulated interaction without free energy */
 641             fscal            = nb_evaluate_single(r2,fr->tab14.scale,fr->tab14.data,qq,c6,c12,&velec,&vvdw);
 642         }
 643
 644         energygrp_elec[gid]  += velec;
 645         energygrp_vdw[gid]   += vvdw;
 646         svmul(fscal,dx,dx);
 647
 648         /* Add the forces */
 649         rvec_inc(f[ai],dx);
 650         rvec_dec(f[aj],dx);
 651
 652         if (g)
 653         {
 654             /* Correct the shift forces using the graph */
 655             ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 656             fshift_index = IVEC2IS(dt);
 657         }
 658         if(fshift_index!=CENTRAL)
 659         {
 660             rvec_inc(fshift[fshift_index],dx);
 661             rvec_dec(fshift[CENTRAL],dx);
 662         }
 663     }
 664     return 0.0;
 665 }
 666
 667