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