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