src/gmxlib/nonbonded/nonbonded.c

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