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 #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, pick a generic one instead */
 289         if (nl->kernelptr_vf == NULL)
 290         {
 291             nl->kernelptr_vf       = (void *) gmx_nb_generic_kernel;
 292             nl->kernelptr_f        = (void *) gmx_nb_generic_kernel;
 293             nl->simd_padding_width = 1;
 294             if (debug)
 295             {
 296                 fprintf(debug,
 297                         "WARNING - Slow generic NB kernel used for neighborlist with\n"
 298                         "    Elec: '%s', Modifier: '%s'\n"
 299                         "    Vdw:  '%s', Modifier: '%s'\n"
 300                         "    Geom: '%s', Other: '%s'\n\n",
 301                         elec, elec_mod, vdw, vdw_mod, geom, other);
 302             }
 303         }
 304     }
 305
 306     return;
 307 }
 308
 309 void do_nonbonded(t_forcerec *fr,
 310                   rvec x[], rvec f_shortrange[], rvec f_longrange[], t_mdatoms *mdatoms, t_blocka *excl,
 311                   gmx_grppairener_t *grppener,
 312                   t_nrnb *nrnb, real *lambda, real *dvdl,
 313                   int nls, int eNL, int flags)
 314 {
 315     t_nblist *        nlist;
 316     int               n, n0, n1, i, i0, i1, sz, range;
 317     t_nblists *       nblists;
 318     nb_kernel_data_t  kernel_data;
 319     nb_kernel_t *     kernelptr = NULL;
 320     rvec *            f;
 321
 322     kernel_data.flags                   = flags;
 323     kernel_data.exclusions              = excl;
 324     kernel_data.lambda                  = lambda;
 325     kernel_data.dvdl                    = dvdl;
 326
 327     if (fr->bAllvsAll)
 328     {
 329         gmx_incons("All-vs-all kernels have not been implemented in version 4.6");
 330         return;
 331     }
 332
 333     if (eNL >= 0)
 334     {
 335         i0 = eNL;
 336         i1 = i0+1;
 337     }
 338     else
 339     {
 340         i0 = 0;
 341         i1 = eNL_NR;
 342     }
 343
 344     if (nls >= 0)
 345     {
 346         n0 = nls;
 347         n1 = nls+1;
 348     }
 349     else
 350     {
 351         n0 = 0;
 352         n1 = fr->nnblists;
 353     }
 354
 355     for (n = n0; (n < n1); n++)
 356     {
 357         nblists = &fr->nblists[n];
 358
 359         kernel_data.table_elec              = &nblists->table_elec;
 360         kernel_data.table_vdw               = &nblists->table_vdw;
 361         kernel_data.table_elec_vdw          = &nblists->table_elec_vdw;
 362
 363         for (range = 0; range < 2; range++)
 364         {
 365             /* Are we doing short/long-range? */
 366             if (range == 0)
 367             {
 368                 /* Short-range */
 369                 if (!(flags & GMX_NONBONDED_DO_SR))
 370                 {
 371                     continue;
 372                 }
 373                 kernel_data.energygrp_elec          = grppener->ener[egCOULSR];
 374                 kernel_data.energygrp_vdw           = grppener->ener[fr->bBHAM ? egBHAMSR : egLJSR];
 375                 kernel_data.energygrp_polarization  = grppener->ener[egGB];
 376                 nlist = nblists->nlist_sr;
 377                 f                                   = f_shortrange;
 378             }
 379             else if (range == 1)
 380             {
 381                 /* Long-range */
 382                 if (!(flags & GMX_NONBONDED_DO_LR))
 383                 {
 384                     continue;
 385                 }
 386                 kernel_data.energygrp_elec          = grppener->ener[egCOULLR];
 387                 kernel_data.energygrp_vdw           = grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR];
 388                 kernel_data.energygrp_polarization  = grppener->ener[egGB];
 389                 nlist = nblists->nlist_lr;
 390                 f                                   = f_longrange;
 391             }
 392
 393             for (i = i0; (i < i1); i++)
 394             {
 395                 if (nlist[i].nri > 0)
 396                 {
 397                     if (flags & GMX_NONBONDED_DO_POTENTIAL)
 398                     {
 399                         /* Potential and force */
 400                         kernelptr = (nb_kernel_t *)nlist[i].kernelptr_vf;
 401                     }
 402                     else
 403                     {
 404                         /* Force only, no potential */
 405                         kernelptr = (nb_kernel_t *)nlist[i].kernelptr_f;
 406                     }
 407
 408                     if (nlist[i].type != GMX_NBLIST_INTERACTION_FREE_ENERGY && (flags & GMX_NONBONDED_DO_FOREIGNLAMBDA))
 409                     {
 410                         /* We don't need the non-perturbed interactions */
 411                         continue;
 412                     }
 413                     (*kernelptr)(&(nlist[i]), x, f, fr, mdatoms, &kernel_data, nrnb);
 414                 }
 415             }
 416         }
 417     }
 418 }
 419
 420 static void
 421 nb_listed_warning_rlimit(const rvec *x, int ai, int aj, int * global_atom_index, real r, real rlimit)
 422 {
 423     gmx_warning("Listed nonbonded interaction between particles %d and %d\n"
 424                 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
 425                 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
 426                 "a smaller molecule you are decoupling during a free energy calculation.\n"
 427                 "Since interactions at distances beyond the table cannot be computed,\n"
 428                 "they are skipped until they are inside the table limit again. You will\n"
 429                 "only see this message once, even if it occurs for several interactions.\n\n"
 430                 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
 431                 "probably something wrong with your system. Only change the table-extension\n"
 432                 "distance in the mdp file if you are really sure that is the reason.\n",
 433                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
 434
 435     if (debug)
 436     {
 437         fprintf(debug,
 438                 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
 439                 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
 440                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
 441     }
 442 }
 443
 444
 445
 446 /* This might logically belong better in the nb_generic.c module, but it is only
 447  * used in do_nonbonded_listed(), and we want it to be inlined there to avoid an
 448  * extra functional call for every single pair listed in the topology.
 449  */
 450 static real
 451 nb_evaluate_single(real r2, real tabscale, real *vftab,
 452                    real qq, real c6, real c12, real *velec, real *vvdw)
 453 {
 454     real       rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
 455     int        ntab;
 456
 457     /* Do the tabulated interactions - first table lookup */
 458     rinv             = gmx_invsqrt(r2);
 459     r                = r2*rinv;
 460     rtab             = r*tabscale;
 461     ntab             = rtab;
 462     eps              = rtab-ntab;
 463     eps2             = eps*eps;
 464     ntab             = 12*ntab;
 465     /* Electrostatics */
 466     Y                = vftab[ntab];
 467     F                = vftab[ntab+1];
 468     Geps             = eps*vftab[ntab+2];
 469     Heps2            = eps2*vftab[ntab+3];
 470     Fp               = F+Geps+Heps2;
 471     VVe              = Y+eps*Fp;
 472     FFe              = Fp+Geps+2.0*Heps2;
 473     /* Dispersion */
 474     Y                = vftab[ntab+4];
 475     F                = vftab[ntab+5];
 476     Geps             = eps*vftab[ntab+6];
 477     Heps2            = eps2*vftab[ntab+7];
 478     Fp               = F+Geps+Heps2;
 479     VVd              = Y+eps*Fp;
 480     FFd              = Fp+Geps+2.0*Heps2;
 481     /* Repulsion */
 482     Y                = vftab[ntab+8];
 483     F                = vftab[ntab+9];
 484     Geps             = eps*vftab[ntab+10];
 485     Heps2            = eps2*vftab[ntab+11];
 486     Fp               = F+Geps+Heps2;
 487     VVr              = Y+eps*Fp;
 488     FFr              = Fp+Geps+2.0*Heps2;
 489
 490     *velec           = qq*VVe;
 491     *vvdw            = c6*VVd+c12*VVr;
 492
 493     fscal            = -(qq*FFe+c6*FFd+c12*FFr)*tabscale*rinv;
 494
 495     return fscal;
 496 }
 497
 498
 499 real
 500 do_nonbonded_listed(int ftype, int nbonds,
 501                     const t_iatom iatoms[], const t_iparams iparams[],
 502                     const rvec x[], rvec f[], rvec fshift[],
 503                     const t_pbc *pbc, const t_graph *g,
 504                     real *lambda, real *dvdl,
 505                     const t_mdatoms *md,
 506                     const t_forcerec *fr, gmx_grppairener_t *grppener,
 507                     int *global_atom_index)
 508 {
 509     int              ielec, ivdw;
 510     real             qq, c6, c12;
 511     rvec             dx;
 512     ivec             dt;
 513     int              i, j, itype, ai, aj, gid;
 514     int              fshift_index;
 515     real             r2, rinv;
 516     real             fscal, velec, vvdw;
 517     real *           energygrp_elec;
 518     real *           energygrp_vdw;
 519     static gmx_bool  warned_rlimit = FALSE;
 520     /* Free energy stuff */
 521     gmx_bool         bFreeEnergy;
 522     real             LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
 523     real             qqB, c6B, c12B, sigma2_def, sigma2_min;
 524
 525
 526     switch (ftype)
 527     {
 528         case F_LJ14:
 529         case F_LJC14_Q:
 530             energygrp_elec = grppener->ener[egCOUL14];
 531             energygrp_vdw  = grppener->ener[egLJ14];
 532             break;
 533         case F_LJC_PAIRS_NB:
 534             energygrp_elec = grppener->ener[egCOULSR];
 535             energygrp_vdw  = grppener->ener[egLJSR];
 536             break;
 537         default:
 538             energygrp_elec = NULL; /* Keep compiler happy */
 539             energygrp_vdw  = NULL; /* Keep compiler happy */
 540             gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
 541             break;
 542     }
 543
 544     if (fr->efep != efepNO)
 545     {
 546         /* Lambda factor for state A=1-lambda and B=lambda */
 547         LFC[0] = 1.0 - lambda[efptCOUL];
 548         LFV[0] = 1.0 - lambda[efptVDW];
 549         LFC[1] = lambda[efptCOUL];
 550         LFV[1] = lambda[efptVDW];
 551
 552         /*derivative of the lambda factor for state A and B */
 553         DLF[0] = -1;
 554         DLF[1] = 1;
 555
 556         /* precalculate */
 557         sigma2_def = pow(fr->sc_sigma6_def, 1.0/3.0);
 558         sigma2_min = pow(fr->sc_sigma6_min, 1.0/3.0);
 559
 560         for (i = 0; i < 2; i++)
 561         {
 562             lfac_coul[i]  = (fr->sc_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 563             dlfac_coul[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFC[i]) : 1);
 564             lfac_vdw[i]   = (fr->sc_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 565             dlfac_vdw[i]  = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFV[i]) : 1);
 566         }
 567     }
 568     else
 569     {
 570         sigma2_min = sigma2_def = 0;
 571     }
 572
 573     bFreeEnergy = FALSE;
 574     for (i = 0; (i < nbonds); )
 575     {
 576         itype = iatoms[i++];
 577         ai    = iatoms[i++];
 578         aj    = iatoms[i++];
 579         gid   = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
 580
 581         /* Get parameters */
 582         switch (ftype)
 583         {
 584             case F_LJ14:
 585                 bFreeEnergy =
 586                     (fr->efep != efepNO &&
 587                      ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
 588                       iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
 589                       iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 590                 qq               = md->chargeA[ai]*md->chargeA[aj]*fr->epsfac*fr->fudgeQQ;
 591                 c6               = iparams[itype].lj14.c6A;
 592                 c12              = iparams[itype].lj14.c12A;
 593                 break;
 594             case F_LJC14_Q:
 595                 qq               = iparams[itype].ljc14.qi*iparams[itype].ljc14.qj*fr->epsfac*iparams[itype].ljc14.fqq;
 596                 c6               = iparams[itype].ljc14.c6;
 597                 c12              = iparams[itype].ljc14.c12;
 598                 break;
 599             case F_LJC_PAIRS_NB:
 600                 qq               = iparams[itype].ljcnb.qi*iparams[itype].ljcnb.qj*fr->epsfac;
 601                 c6               = iparams[itype].ljcnb.c6;
 602                 c12              = iparams[itype].ljcnb.c12;
 603                 break;
 604             default:
 605                 /* Cannot happen since we called gmx_fatal() above in this case */
 606                 qq = c6 = c12 = 0; /* Keep compiler happy */
 607                 break;
 608         }
 609
 610         /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
 611          * included in the general nfbp array now. This means the tables are scaled down by the
 612          * same factor, so when we use the original c6/c12 parameters from iparams[] they must
 613          * be scaled up.
 614          */
 615         c6  *= 6.0;
 616         c12 *= 12.0;
 617
 618         /* Do we need to apply full periodic boundary conditions? */
 619         if (fr->bMolPBC == TRUE)
 620         {
 621             fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 622         }
 623         else
 624         {
 625             fshift_index = CENTRAL;
 626             rvec_sub(x[ai], x[aj], dx);
 627         }
 628         r2           = norm2(dx);
 629
 630         if (r2 >= fr->tab14.r*fr->tab14.r)
 631         {
 632             if (warned_rlimit == FALSE)
 633             {
 634                 nb_listed_warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->tab14.r);
 635                 warned_rlimit = TRUE;
 636             }
 637             continue;
 638         }
 639
 640         if (bFreeEnergy)
 641         {
 642             /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
 643             qqB              = md->chargeB[ai]*md->chargeB[aj]*fr->epsfac*fr->fudgeQQ;
 644             c6B              = iparams[itype].lj14.c6B*6.0;
 645             c12B             = iparams[itype].lj14.c12B*12.0;
 646
 647             fscal            = nb_free_energy_evaluate_single(r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw,
 648                                                               fr->tab14.scale, fr->tab14.data, qq, c6, c12, qqB, c6B, c12B,
 649                                                               LFC, LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw,
 650                                                               fr->sc_sigma6_def, fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
 651         }
 652         else
 653         {
 654             /* Evaluate tabulated interaction without free energy */
 655             fscal            = nb_evaluate_single(r2, fr->tab14.scale, fr->tab14.data, qq, c6, c12, &velec, &vvdw);
 656         }
 657
 658         energygrp_elec[gid]  += velec;
 659         energygrp_vdw[gid]   += vvdw;
 660         svmul(fscal, dx, dx);
 661
 662         /* Add the forces */
 663         rvec_inc(f[ai], dx);
 664         rvec_dec(f[aj], dx);
 665
 666         if (g)
 667         {
 668             /* Correct the shift forces using the graph */
 669             ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
 670             fshift_index = IVEC2IS(dt);
 671         }
 672         if (fshift_index != CENTRAL)
 673         {
 674             rvec_inc(fshift[fshift_index], dx);
 675             rvec_dec(fshift[CENTRAL], dx);
 676         }
 677     }
 678     return 0.0;
 679 }