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