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