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