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