src/gmxlib/nonbonded/nonbonded.c

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