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