src/gromacs/gmxlib/nonbonded/nb_free_energy.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 #include "gmxpre.h"
  38
  39 #include "nb_free_energy.h"
  40
  41 #include <math.h>
  42
  43 #include "gromacs/gmxlib/nonbonded/nb_kernel.h"
  44 #include "gromacs/legacyheaders/macros.h"
  45 #include "gromacs/legacyheaders/nonbonded.h"
  46 #include "gromacs/legacyheaders/nrnb.h"
  47 #include "gromacs/legacyheaders/typedefs.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/utility/fatalerror.h"
  50
  51 void
  52 gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
  53                           rvec * gmx_restrict              xx,
  54                           rvec * gmx_restrict              ff,
  55                           t_forcerec * gmx_restrict        fr,
  56                           const t_mdatoms * gmx_restrict   mdatoms,
  57                           nb_kernel_data_t * gmx_restrict  kernel_data,
  58                           t_nrnb * gmx_restrict            nrnb)
  59 {
  60
  61 #define  STATE_A  0
  62 #define  STATE_B  1
  63 #define  NSTATES  2
  64     int           i, j, n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid;
  65     real          shX, shY, shZ;
  66     real          tx, ty, tz, Fscal;
  67     double        FscalC[NSTATES], FscalV[NSTATES];  /* Needs double for sc_power==48 */
  68     double        Vcoul[NSTATES], Vvdw[NSTATES];     /* Needs double for sc_power==48 */
  69     real          rinv6, r, rt, rtC, rtV;
  70     real          iqA, iqB;
  71     real          qq[NSTATES], vctot, krsq;
  72     int           ntiA, ntiB, tj[NSTATES];
  73     real          Vvdw6, Vvdw12, vvtot;
  74     real          ix, iy, iz, fix, fiy, fiz;
  75     real          dx, dy, dz, rsq, rinv;
  76     real          c6[NSTATES], c12[NSTATES], c6grid;
  77     real          LFC[NSTATES], LFV[NSTATES], DLF[NSTATES];
  78     double        dvdl_coul, dvdl_vdw;
  79     real          lfac_coul[NSTATES], dlfac_coul[NSTATES], lfac_vdw[NSTATES], dlfac_vdw[NSTATES];
  80     real          sigma6[NSTATES], alpha_vdw_eff, alpha_coul_eff, sigma2_def, sigma2_min;
  81     double        rp, rpm2, rC, rV, rinvC, rpinvC, rinvV, rpinvV; /* Needs double for sc_power==48 */
  82     real          sigma2[NSTATES], sigma_pow[NSTATES], sigma_powm2[NSTATES], rs, rs2;
  83     int           do_tab, tab_elemsize;
  84     int           n0, n1C, n1V, nnn;
  85     real          Y, F, G, H, Fp, Geps, Heps2, epsC, eps2C, epsV, eps2V, VV, FF;
  86     int           icoul, ivdw;
  87     int           nri;
  88     const int *   iinr;
  89     const int *   jindex;
  90     const int *   jjnr;
  91     const int *   shift;
  92     const int *   gid;
  93     const int *   typeA;
  94     const int *   typeB;
  95     int           ntype;
  96     const real *  shiftvec;
  97     real          dvdl_part;
  98     real *        fshift;
  99     real          tabscale = 0;
 100     const real *  VFtab    = NULL;
 101     const real *  x;
 102     real *        f;
 103     real          facel, krf, crf;
 104     const real *  chargeA;
 105     const real *  chargeB;
 106     real          sigma6_min, sigma6_def, lam_power, sc_power, sc_r_power;
 107     real          alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc_lj;
 108     real          ewcljrsq, ewclj, ewclj2, exponent, poly, vvdw_disp, vvdw_rep, sh_lj_ewald;
 109     real          ewclj6;
 110     const real *  nbfp, *nbfp_grid;
 111     real *        dvdl;
 112     real *        Vv;
 113     real *        Vc;
 114     gmx_bool      bDoForces, bDoShiftForces, bDoPotential;
 115     real          rcoulomb, rvdw, sh_invrc6;
 116     gmx_bool      bExactElecCutoff, bExactVdwCutoff, bExactCutoffAll;
 117     gmx_bool      bEwald, bEwaldLJ;
 118     real          rcutoff_max2;
 119     const real *  tab_ewald_F_lj;
 120     const real *  tab_ewald_V_lj;
 121     real          d, d2, sw, dsw, rinvcorr;
 122     real          elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
 123     real          vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
 124     gmx_bool      bConvertEwaldToCoulomb, bConvertLJEwaldToLJ6;
 125     gmx_bool      bComputeVdwInteraction, bComputeElecInteraction;
 126     const real *  ewtab;
 127     int           ewitab;
 128     real          ewrt, eweps, ewtabscale, ewtabhalfspace, sh_ewald;
 129
 130     const real    onetwelfth  = 1.0/12.0;
 131     const real    onesixth    = 1.0/6.0;
 132     const real    zero        = 0.0;
 133     const real    half        = 0.5;
 134     const real    one         = 1.0;
 135     const real    two         = 2.0;
 136     const real    six         = 6.0;
 137     const real    fourtyeight = 48.0;
 138
 139     sh_ewald            = fr->ic->sh_ewald;
 140     ewtab               = fr->ic->tabq_coul_FDV0;
 141     ewtabscale          = fr->ic->tabq_scale;
 142     ewtabhalfspace      = half/ewtabscale;
 143     tab_ewald_F_lj      = fr->ic->tabq_vdw_F;
 144     tab_ewald_V_lj      = fr->ic->tabq_vdw_V;
 145
 146     x                   = xx[0];
 147     f                   = ff[0];
 148
 149     fshift              = fr->fshift[0];
 150
 151     nri                 = nlist->nri;
 152     iinr                = nlist->iinr;
 153     jindex              = nlist->jindex;
 154     jjnr                = nlist->jjnr;
 155     icoul               = nlist->ielec;
 156     ivdw                = nlist->ivdw;
 157     shift               = nlist->shift;
 158     gid                 = nlist->gid;
 159
 160     shiftvec            = fr->shift_vec[0];
 161     chargeA             = mdatoms->chargeA;
 162     chargeB             = mdatoms->chargeB;
 163     facel               = fr->epsfac;
 164     krf                 = fr->k_rf;
 165     crf                 = fr->c_rf;
 166     ewc_lj              = fr->ewaldcoeff_lj;
 167     Vc                  = kernel_data->energygrp_elec;
 168     typeA               = mdatoms->typeA;
 169     typeB               = mdatoms->typeB;
 170     ntype               = fr->ntype;
 171     nbfp                = fr->nbfp;
 172     nbfp_grid           = fr->ljpme_c6grid;
 173     Vv                  = kernel_data->energygrp_vdw;
 174     lambda_coul         = kernel_data->lambda[efptCOUL];
 175     lambda_vdw          = kernel_data->lambda[efptVDW];
 176     dvdl                = kernel_data->dvdl;
 177     alpha_coul          = fr->sc_alphacoul;
 178     alpha_vdw           = fr->sc_alphavdw;
 179     lam_power           = fr->sc_power;
 180     sc_r_power          = fr->sc_r_power;
 181     sigma6_def          = fr->sc_sigma6_def;
 182     sigma6_min          = fr->sc_sigma6_min;
 183     bDoForces           = kernel_data->flags & GMX_NONBONDED_DO_FORCE;
 184     bDoShiftForces      = kernel_data->flags & GMX_NONBONDED_DO_SHIFTFORCE;
 185     bDoPotential        = kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL;
 186
 187     rcoulomb            = fr->rcoulomb;
 188     rvdw                = fr->rvdw;
 189     sh_invrc6           = fr->ic->sh_invrc6;
 190     sh_lj_ewald         = fr->ic->sh_lj_ewald;
 191     ewclj               = fr->ewaldcoeff_lj;
 192     ewclj2              = ewclj*ewclj;
 193     ewclj6              = ewclj2*ewclj2*ewclj2;
 194
 195     if (fr->coulomb_modifier == eintmodPOTSWITCH)
 196     {
 197         d               = fr->rcoulomb-fr->rcoulomb_switch;
 198         elec_swV3       = -10.0/(d*d*d);
 199         elec_swV4       =  15.0/(d*d*d*d);
 200         elec_swV5       =  -6.0/(d*d*d*d*d);
 201         elec_swF2       = -30.0/(d*d*d);
 202         elec_swF3       =  60.0/(d*d*d*d);
 203         elec_swF4       = -30.0/(d*d*d*d*d);
 204     }
 205     else
 206     {
 207         /* Avoid warnings from stupid compilers (looking at you, Clang!) */
 208         elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
 209     }
 210
 211     if (fr->vdw_modifier == eintmodPOTSWITCH)
 212     {
 213         d               = fr->rvdw-fr->rvdw_switch;
 214         vdw_swV3        = -10.0/(d*d*d);
 215         vdw_swV4        =  15.0/(d*d*d*d);
 216         vdw_swV5        =  -6.0/(d*d*d*d*d);
 217         vdw_swF2        = -30.0/(d*d*d);
 218         vdw_swF3        =  60.0/(d*d*d*d);
 219         vdw_swF4        = -30.0/(d*d*d*d*d);
 220     }
 221     else
 222     {
 223         /* Avoid warnings from stupid compilers (looking at you, Clang!) */
 224         vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
 225     }
 226
 227     if (fr->cutoff_scheme == ecutsVERLET)
 228     {
 229         const interaction_const_t *ic;
 230
 231         ic = fr->ic;
 232         if (EVDW_PME(ic->vdwtype))
 233         {
 234             ivdw         = GMX_NBKERNEL_VDW_LJEWALD;
 235         }
 236         else
 237         {
 238             ivdw         = GMX_NBKERNEL_VDW_LENNARDJONES;
 239         }
 240
 241         if (ic->eeltype == eelCUT || EEL_RF(ic->eeltype))
 242         {
 243             icoul        = GMX_NBKERNEL_ELEC_REACTIONFIELD;
 244         }
 245         else if (EEL_PME_EWALD(ic->eeltype))
 246         {
 247             icoul        = GMX_NBKERNEL_ELEC_EWALD;
 248         }
 249         else
 250         {
 251             gmx_incons("Unsupported eeltype with Verlet and free-energy");
 252         }
 253
 254         bExactElecCutoff = TRUE;
 255         bExactVdwCutoff  = TRUE;
 256     }
 257     else
 258     {
 259         bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
 260         bExactVdwCutoff  = (fr->vdw_modifier != eintmodNONE);
 261     }
 262
 263     bExactCutoffAll = (bExactElecCutoff && bExactVdwCutoff);
 264     rcutoff_max2    = max(fr->rcoulomb, fr->rvdw);
 265     rcutoff_max2    = rcutoff_max2*rcutoff_max2;
 266
 267     bEwald          = (icoul == GMX_NBKERNEL_ELEC_EWALD);
 268     bEwaldLJ        = (ivdw == GMX_NBKERNEL_VDW_LJEWALD);
 269
 270     /* For Ewald/PME interactions we cannot easily apply the soft-core component to
 271      * reciprocal space. When we use vanilla (not switch/shift) Ewald interactions, we
 272      * can apply the small trick of subtracting the _reciprocal_ space contribution
 273      * in this kernel, and instead apply the free energy interaction to the 1/r
 274      * (standard coulomb) interaction.
 275      *
 276      * However, we cannot use this approach for switch-modified since we would then
 277      * effectively end up evaluating a significantly different interaction here compared to the
 278      * normal (non-free-energy) kernels, either by applying a cutoff at a different
 279      * position than what the user requested, or by switching different
 280      * things (1/r rather than short-range Ewald). For these settings, we just
 281      * use the traditional short-range Ewald interaction in that case.
 282      */
 283     bConvertEwaldToCoulomb = (bEwald && (fr->coulomb_modifier != eintmodPOTSWITCH));
 284     /* For now the below will always be true (since LJ-PME only works with Shift in Gromacs-5.0),
 285      * but writing it this way means we stay in sync with coulomb, and it avoids future bugs.
 286      */
 287     bConvertLJEwaldToLJ6   = (bEwaldLJ && (fr->vdw_modifier   != eintmodPOTSWITCH));
 288
 289     /* We currently don't implement exclusion correction, needed with the Verlet cut-off scheme, without conversion */
 290     if (fr->cutoff_scheme == ecutsVERLET &&
 291         ((bEwald   && !bConvertEwaldToCoulomb) ||
 292          (bEwaldLJ && !bConvertLJEwaldToLJ6)))
 293     {
 294         gmx_incons("Unimplemented non-bonded setup");
 295     }
 296
 297     /* fix compiler warnings */
 298     nj1   = 0;
 299     n1C   = n1V   = 0;
 300     epsC  = epsV  = 0;
 301     eps2C = eps2V = 0;
 302
 303     dvdl_coul  = 0;
 304     dvdl_vdw   = 0;
 305
 306     /* Lambda factor for state A, 1-lambda*/
 307     LFC[STATE_A] = one - lambda_coul;
 308     LFV[STATE_A] = one - lambda_vdw;
 309
 310     /* Lambda factor for state B, lambda*/
 311     LFC[STATE_B] = lambda_coul;
 312     LFV[STATE_B] = lambda_vdw;
 313
 314     /*derivative of the lambda factor for state A and B */
 315     DLF[STATE_A] = -1;
 316     DLF[STATE_B] = 1;
 317
 318     for (i = 0; i < NSTATES; i++)
 319     {
 320         lfac_coul[i]  = (lam_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 321         dlfac_coul[i] = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFC[i]) : 1);
 322         lfac_vdw[i]   = (lam_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 323         dlfac_vdw[i]  = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFV[i]) : 1);
 324     }
 325     /* precalculate */
 326     sigma2_def = pow(sigma6_def, 1.0/3.0);
 327     sigma2_min = pow(sigma6_min, 1.0/3.0);
 328
 329     /* Ewald (not PME) table is special (icoul==enbcoulFEWALD) */
 330
 331     do_tab = (icoul == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE ||
 332               ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE);
 333     if (do_tab)
 334     {
 335         tabscale         = kernel_data->table_elec_vdw->scale;
 336         VFtab            = kernel_data->table_elec_vdw->data;
 337         /* we always use the combined table here */
 338         tab_elemsize     = 12;
 339     }
 340
 341     for (n = 0; (n < nri); n++)
 342     {
 343         int npair_within_cutoff;
 344
 345         npair_within_cutoff = 0;
 346
 347         is3              = 3*shift[n];
 348         shX              = shiftvec[is3];
 349         shY              = shiftvec[is3+1];
 350         shZ              = shiftvec[is3+2];
 351         nj0              = jindex[n];
 352         nj1              = jindex[n+1];
 353         ii               = iinr[n];
 354         ii3              = 3*ii;
 355         ix               = shX + x[ii3+0];
 356         iy               = shY + x[ii3+1];
 357         iz               = shZ + x[ii3+2];
 358         iqA              = facel*chargeA[ii];
 359         iqB              = facel*chargeB[ii];
 360         ntiA             = 2*ntype*typeA[ii];
 361         ntiB             = 2*ntype*typeB[ii];
 362         vctot            = 0;
 363         vvtot            = 0;
 364         fix              = 0;
 365         fiy              = 0;
 366         fiz              = 0;
 367
 368         for (k = nj0; (k < nj1); k++)
 369         {
 370             jnr              = jjnr[k];
 371             j3               = 3*jnr;
 372             dx               = ix - x[j3];
 373             dy               = iy - x[j3+1];
 374             dz               = iz - x[j3+2];
 375             rsq              = dx*dx + dy*dy + dz*dz;
 376
 377             if (bExactCutoffAll && rsq >= rcutoff_max2)
 378             {
 379                 /* We save significant time by skipping all code below.
 380                  * Note that with soft-core interactions, the actual cut-off
 381                  * check might be different. But since the soft-core distance
 382                  * is always larger than r, checking on r here is safe.
 383                  */
 384                 continue;
 385             }
 386             npair_within_cutoff++;
 387
 388             if (rsq > 0)
 389             {
 390                 rinv         = gmx_invsqrt(rsq);
 391                 r            = rsq*rinv;
 392             }
 393             else
 394             {
 395                 /* The force at r=0 is zero, because of symmetry.
 396                  * But note that the potential is in general non-zero,
 397                  * since the soft-cored r will be non-zero.
 398                  */
 399                 rinv         = 0;
 400                 r            = 0;
 401             }
 402
 403             if (sc_r_power == six)
 404             {
 405                 rpm2             = rsq*rsq;  /* r4 */
 406                 rp               = rpm2*rsq; /* r6 */
 407             }
 408             else if (sc_r_power == fourtyeight)
 409             {
 410                 rp               = rsq*rsq*rsq; /* r6 */
 411                 rp               = rp*rp;       /* r12 */
 412                 rp               = rp*rp;       /* r24 */
 413                 rp               = rp*rp;       /* r48 */
 414                 rpm2             = rp/rsq;      /* r46 */
 415             }
 416             else
 417             {
 418                 rp             = pow(r, sc_r_power);  /* not currently supported as input, but can handle it */
 419                 rpm2           = rp/rsq;
 420             }
 421
 422             Fscal = 0;
 423
 424             qq[STATE_A]      = iqA*chargeA[jnr];
 425             qq[STATE_B]      = iqB*chargeB[jnr];
 426
 427             tj[STATE_A]      = ntiA+2*typeA[jnr];
 428             tj[STATE_B]      = ntiB+2*typeB[jnr];
 429
 430             if (nlist->excl_fep == NULL || nlist->excl_fep[k])
 431             {
 432                 c6[STATE_A]      = nbfp[tj[STATE_A]];
 433                 c6[STATE_B]      = nbfp[tj[STATE_B]];
 434
 435                 for (i = 0; i < NSTATES; i++)
 436                 {
 437                     c12[i]             = nbfp[tj[i]+1];
 438                     if ((c6[i] > 0) && (c12[i] > 0))
 439                     {
 440                         /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
 441                         sigma6[i]       = half*c12[i]/c6[i];
 442                         sigma2[i]       = pow(sigma6[i], 1.0/3.0);
 443                         /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 444                            what data to store externally.  Can't be fixed without larger scale changes, so not 4.6 */
 445                         if (sigma6[i] < sigma6_min)   /* for disappearing coul and vdw with soft core at the same time */
 446                         {
 447                             sigma6[i] = sigma6_min;
 448                             sigma2[i] = sigma2_min;
 449                         }
 450                     }
 451                     else
 452                     {
 453                         sigma6[i]       = sigma6_def;
 454                         sigma2[i]       = sigma2_def;
 455                     }
 456                     if (sc_r_power == six)
 457                     {
 458                         sigma_pow[i]    = sigma6[i];
 459                         sigma_powm2[i]  = sigma6[i]/sigma2[i];
 460                     }
 461                     else if (sc_r_power == fourtyeight)
 462                     {
 463                         sigma_pow[i]    = sigma6[i]*sigma6[i];       /* sigma^12 */
 464                         sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
 465                         sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
 466                         sigma_powm2[i]  = sigma_pow[i]/sigma2[i];
 467                     }
 468                     else
 469                     {    /* not really supported as input, but in here for testing the general case*/
 470                         sigma_pow[i]    = pow(sigma2[i], sc_r_power/2);
 471                         sigma_powm2[i]  = sigma_pow[i]/(sigma2[i]);
 472                     }
 473                 }
 474
 475                 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 476                 if ((c12[STATE_A] > 0) && (c12[STATE_B] > 0))
 477                 {
 478                     alpha_vdw_eff    = 0;
 479                     alpha_coul_eff   = 0;
 480                 }
 481                 else
 482                 {
 483                     alpha_vdw_eff    = alpha_vdw;
 484                     alpha_coul_eff   = alpha_coul;
 485                 }
 486
 487                 for (i = 0; i < NSTATES; i++)
 488                 {
 489                     FscalC[i]    = 0;
 490                     FscalV[i]    = 0;
 491                     Vcoul[i]     = 0;
 492                     Vvdw[i]      = 0;
 493
 494                     /* Only spend time on A or B state if it is non-zero */
 495                     if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
 496                     {
 497                         /* this section has to be inside the loop because of the dependence on sigma_pow */
 498                         rpinvC         = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
 499                         rinvC          = pow(rpinvC, one/sc_r_power);
 500                         rC             = one/rinvC;
 501
 502                         rpinvV         = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
 503                         rinvV          = pow(rpinvV, one/sc_r_power);
 504                         rV             = one/rinvV;
 505
 506                         if (do_tab)
 507                         {
 508                             rtC        = rC*tabscale;
 509                             n0         = rtC;
 510                             epsC       = rtC-n0;
 511                             eps2C      = epsC*epsC;
 512                             n1C        = tab_elemsize*n0;
 513
 514                             rtV        = rV*tabscale;
 515                             n0         = rtV;
 516                             epsV       = rtV-n0;
 517                             eps2V      = epsV*epsV;
 518                             n1V        = tab_elemsize*n0;
 519                         }
 520
 521                         /* Only process the coulomb interactions if we have charges,
 522                          * and if we either include all entries in the list (no cutoff
 523                          * used in the kernel), or if we are within the cutoff.
 524                          */
 525                         bComputeElecInteraction = !bExactElecCutoff ||
 526                             ( bConvertEwaldToCoulomb && r < rcoulomb) ||
 527                             (!bConvertEwaldToCoulomb && rC < rcoulomb);
 528
 529                         if ( (qq[i] != 0) && bComputeElecInteraction)
 530                         {
 531                             switch (icoul)
 532                             {
 533                                 case GMX_NBKERNEL_ELEC_COULOMB:
 534                                     /* simple cutoff */
 535                                     Vcoul[i]   = qq[i]*rinvC;
 536                                     FscalC[i]  = Vcoul[i];
 537                                     /* The shift for the Coulomb potential is stored in
 538                                      * the RF parameter c_rf, which is 0 without shift.
 539                                      */
 540                                     Vcoul[i]  -= qq[i]*fr->ic->c_rf;
 541                                     break;
 542
 543                                 case GMX_NBKERNEL_ELEC_REACTIONFIELD:
 544                                     /* reaction-field */
 545                                     Vcoul[i]   = qq[i]*(rinvC + krf*rC*rC-crf);
 546                                     FscalC[i]  = qq[i]*(rinvC - two*krf*rC*rC);
 547                                     break;
 548
 549                                 case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE:
 550                                     /* non-Ewald tabulated coulomb */
 551                                     nnn        = n1C;
 552                                     Y          = VFtab[nnn];
 553                                     F          = VFtab[nnn+1];
 554                                     Geps       = epsC*VFtab[nnn+2];
 555                                     Heps2      = eps2C*VFtab[nnn+3];
 556                                     Fp         = F+Geps+Heps2;
 557                                     VV         = Y+epsC*Fp;
 558                                     FF         = Fp+Geps+two*Heps2;
 559                                     Vcoul[i]   = qq[i]*VV;
 560                                     FscalC[i]  = -qq[i]*tabscale*FF*rC;
 561                                     break;
 562
 563                                 case GMX_NBKERNEL_ELEC_GENERALIZEDBORN:
 564                                     gmx_fatal(FARGS, "Free energy and GB not implemented.\n");
 565                                     break;
 566
 567                                 case GMX_NBKERNEL_ELEC_EWALD:
 568                                     if (bConvertEwaldToCoulomb)
 569                                     {
 570                                         /* Ewald FEP is done only on the 1/r part */
 571                                         Vcoul[i]   = qq[i]*(rinvC-sh_ewald);
 572                                         FscalC[i]  = qq[i]*rinvC;
 573                                     }
 574                                     else
 575                                     {
 576                                         ewrt      = rC*ewtabscale;
 577                                         ewitab    = (int) ewrt;
 578                                         eweps     = ewrt-ewitab;
 579                                         ewitab    = 4*ewitab;
 580                                         FscalC[i] = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 581                                         rinvcorr  = rinvC-sh_ewald;
 582                                         Vcoul[i]  = qq[i]*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+FscalC[i])));
 583                                         FscalC[i] = qq[i]*(rinvC-rC*FscalC[i]);
 584                                     }
 585                                     break;
 586
 587                                 case GMX_NBKERNEL_ELEC_NONE:
 588                                     FscalC[i]  = zero;
 589                                     Vcoul[i]   = zero;
 590                                     break;
 591
 592                                 default:
 593                                     gmx_incons("Invalid icoul in free energy kernel");
 594                                     break;
 595                             }
 596
 597                             if (fr->coulomb_modifier == eintmodPOTSWITCH)
 598                             {
 599                                 d                = rC-fr->rcoulomb_switch;
 600                                 d                = (d > zero) ? d : zero;
 601                                 d2               = d*d;
 602                                 sw               = one+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
 603                                 dsw              = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
 604
 605                                 FscalC[i]        = FscalC[i]*sw - rC*Vcoul[i]*dsw;
 606                                 Vcoul[i]        *= sw;
 607
 608                                 FscalC[i]        = (rC < rcoulomb) ? FscalC[i] : zero;
 609                                 Vcoul[i]         = (rC < rcoulomb) ? Vcoul[i] : zero;
 610                             }
 611                         }
 612
 613                         /* Only process the VDW interactions if we have
 614                          * some non-zero parameters, and if we either
 615                          * include all entries in the list (no cutoff used
 616                          * in the kernel), or if we are within the cutoff.
 617                          */
 618                         bComputeVdwInteraction = !bExactVdwCutoff ||
 619                             ( bConvertLJEwaldToLJ6 && r < rvdw) ||
 620                             (!bConvertLJEwaldToLJ6 && rV < rvdw);
 621                         if ((c6[i] != 0 || c12[i] != 0) && bComputeVdwInteraction)
 622                         {
 623                             switch (ivdw)
 624                             {
 625                                 case GMX_NBKERNEL_VDW_LENNARDJONES:
 626                                     /* cutoff LJ */
 627                                     if (sc_r_power == six)
 628                                     {
 629                                         rinv6            = rpinvV;
 630                                     }
 631                                     else
 632                                     {
 633                                         rinv6            = rinvV*rinvV;
 634                                         rinv6            = rinv6*rinv6*rinv6;
 635                                     }
 636                                     Vvdw6            = c6[i]*rinv6;
 637                                     Vvdw12           = c12[i]*rinv6*rinv6;
 638
 639                                     Vvdw[i]          = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*onetwelfth
 640                                                          - (Vvdw6 - c6[i]*sh_invrc6)*onesixth);
 641                                     FscalV[i]        = Vvdw12 - Vvdw6;
 642                                     break;
 643
 644                                 case GMX_NBKERNEL_VDW_BUCKINGHAM:
 645                                     gmx_fatal(FARGS, "Buckingham free energy not supported.");
 646                                     break;
 647
 648                                 case GMX_NBKERNEL_VDW_CUBICSPLINETABLE:
 649                                     /* Table LJ */
 650                                     nnn = n1V+4;
 651                                     /* dispersion */
 652                                     Y          = VFtab[nnn];
 653                                     F          = VFtab[nnn+1];
 654                                     Geps       = epsV*VFtab[nnn+2];
 655                                     Heps2      = eps2V*VFtab[nnn+3];
 656                                     Fp         = F+Geps+Heps2;
 657                                     VV         = Y+epsV*Fp;
 658                                     FF         = Fp+Geps+two*Heps2;
 659                                     Vvdw[i]   += c6[i]*VV;
 660                                     FscalV[i] -= c6[i]*tabscale*FF*rV;
 661
 662                                     /* repulsion */
 663                                     Y          = VFtab[nnn+4];
 664                                     F          = VFtab[nnn+5];
 665                                     Geps       = epsV*VFtab[nnn+6];
 666                                     Heps2      = eps2V*VFtab[nnn+7];
 667                                     Fp         = F+Geps+Heps2;
 668                                     VV         = Y+epsV*Fp;
 669                                     FF         = Fp+Geps+two*Heps2;
 670                                     Vvdw[i]   += c12[i]*VV;
 671                                     FscalV[i] -= c12[i]*tabscale*FF*rV;
 672                                     break;
 673
 674                                 case GMX_NBKERNEL_VDW_LJEWALD:
 675                                     if (sc_r_power == six)
 676                                     {
 677                                         rinv6            = rpinvV;
 678                                     }
 679                                     else
 680                                     {
 681                                         rinv6            = rinvV*rinvV;
 682                                         rinv6            = rinv6*rinv6*rinv6;
 683                                     }
 684                                     c6grid           = nbfp_grid[tj[i]];
 685
 686                                     if (bConvertLJEwaldToLJ6)
 687                                     {
 688                                         /* cutoff LJ */
 689                                         Vvdw6            = c6[i]*rinv6;
 690                                         Vvdw12           = c12[i]*rinv6*rinv6;
 691
 692                                         Vvdw[i]          = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*onetwelfth
 693                                                              - (Vvdw6 - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)*onesixth);
 694                                         FscalV[i]        = Vvdw12 - Vvdw6;
 695                                     }
 696                                     else
 697                                     {
 698                                         /* Normal LJ-PME */
 699                                         ewcljrsq         = ewclj2*rV*rV;
 700                                         exponent         = exp(-ewcljrsq);
 701                                         poly             = exponent*(one + ewcljrsq + ewcljrsq*ewcljrsq*half);
 702                                         vvdw_disp        = (c6[i]-c6grid*(one-poly))*rinv6;
 703                                         vvdw_rep         = c12[i]*rinv6*rinv6;
 704                                         FscalV[i]        = vvdw_rep - vvdw_disp - c6grid*onesixth*exponent*ewclj6;
 705                                         Vvdw[i]          = (vvdw_rep - c12[i]*sh_invrc6*sh_invrc6)*onetwelfth - (vvdw_disp - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)/six;
 706                                     }
 707                                     break;
 708
 709                                 case GMX_NBKERNEL_VDW_NONE:
 710                                     Vvdw[i]    = zero;
 711                                     FscalV[i]  = zero;
 712                                     break;
 713
 714                                 default:
 715                                     gmx_incons("Invalid ivdw in free energy kernel");
 716                                     break;
 717                             }
 718
 719                             if (fr->vdw_modifier == eintmodPOTSWITCH)
 720                             {
 721                                 d                = rV-fr->rvdw_switch;
 722                                 d                = (d > zero) ? d : zero;
 723                                 d2               = d*d;
 724                                 sw               = one+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
 725                                 dsw              = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
 726
 727                                 FscalV[i]        = FscalV[i]*sw - rV*Vvdw[i]*dsw;
 728                                 Vvdw[i]         *= sw;
 729
 730                                 FscalV[i]  = (rV < rvdw) ? FscalV[i] : zero;
 731                                 Vvdw[i]    = (rV < rvdw) ? Vvdw[i] : zero;
 732                             }
 733                         }
 734
 735                         /* FscalC (and FscalV) now contain: dV/drC * rC
 736                          * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p
 737                          * Further down we first multiply by r^p-2 and then by
 738                          * the vector r, which in total gives: dV/drC * (r/rC)^1-p
 739                          */
 740                         FscalC[i] *= rpinvC;
 741                         FscalV[i] *= rpinvV;
 742                     }
 743                 }
 744
 745                 /* Assemble A and B states */
 746                 for (i = 0; i < NSTATES; i++)
 747                 {
 748                     vctot         += LFC[i]*Vcoul[i];
 749                     vvtot         += LFV[i]*Vvdw[i];
 750
 751                     Fscal         += LFC[i]*FscalC[i]*rpm2;
 752                     Fscal         += LFV[i]*FscalV[i]*rpm2;
 753
 754                     dvdl_coul     += Vcoul[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*FscalC[i]*sigma_pow[i];
 755                     dvdl_vdw      += Vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*FscalV[i]*sigma_pow[i];
 756                 }
 757             }
 758             else if (icoul == GMX_NBKERNEL_ELEC_REACTIONFIELD)
 759             {
 760                 /* For excluded pairs, which are only in this pair list when
 761                  * using the Verlet scheme, we don't use soft-core.
 762                  * The group scheme also doesn't soft-core for these.
 763                  * As there is no singularity, there is no need for soft-core.
 764                  */
 765                 VV = krf*rsq - crf;
 766                 FF = -two*krf;
 767
 768                 if (ii == jnr)
 769                 {
 770                     VV *= half;
 771                 }
 772
 773                 for (i = 0; i < NSTATES; i++)
 774                 {
 775                     vctot      += LFC[i]*qq[i]*VV;
 776                     Fscal      += LFC[i]*qq[i]*FF;
 777                     dvdl_coul  += DLF[i]*qq[i]*VV;
 778                 }
 779             }
 780
 781             if (bConvertEwaldToCoulomb && ( !bExactElecCutoff || r < rcoulomb ) )
 782             {
 783                 /* See comment in the preamble. When using Ewald interactions
 784                  * (unless we use a switch modifier) we subtract the reciprocal-space
 785                  * Ewald component here which made it possible to apply the free
 786                  * energy interaction to 1/r (vanilla coulomb short-range part)
 787                  * above. This gets us closer to the ideal case of applying
 788                  * the softcore to the entire electrostatic interaction,
 789                  * including the reciprocal-space component.
 790                  */
 791                 real v_lr, f_lr;
 792
 793                 ewrt      = r*ewtabscale;
 794                 ewitab    = (int) ewrt;
 795                 eweps     = ewrt-ewitab;
 796                 ewitab    = 4*ewitab;
 797                 f_lr      = ewtab[ewitab]+eweps*ewtab[ewitab+1];
 798                 v_lr      = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr));
 799                 f_lr     *= rinv;
 800
 801                 /* Note that any possible Ewald shift has already been applied in
 802                  * the normal interaction part above.
 803                  */
 804
 805                 if (ii == jnr)
 806                 {
 807                     /* If we get here, the i particle (ii) has itself (jnr)
 808                      * in its neighborlist. This can only happen with the Verlet
 809                      * scheme, and corresponds to a self-interaction that will
 810                      * occur twice. Scale it down by 50% to only include it once.
 811                      */
 812                     v_lr *= half;
 813                 }
 814
 815                 for (i = 0; i < NSTATES; i++)
 816                 {
 817                     vctot      -= LFC[i]*qq[i]*v_lr;
 818                     Fscal      -= LFC[i]*qq[i]*f_lr;
 819                     dvdl_coul  -= (DLF[i]*qq[i])*v_lr;
 820                 }
 821             }
 822
 823             if (bConvertLJEwaldToLJ6 && (!bExactVdwCutoff || r < rvdw))
 824             {
 825                 /* See comment in the preamble. When using LJ-Ewald interactions
 826                  * (unless we use a switch modifier) we subtract the reciprocal-space
 827                  * Ewald component here which made it possible to apply the free
 828                  * energy interaction to r^-6 (vanilla LJ6 short-range part)
 829                  * above. This gets us closer to the ideal case of applying
 830                  * the softcore to the entire VdW interaction,
 831                  * including the reciprocal-space component.
 832                  */
 833                 /* We could also use the analytical form here
 834                  * iso a table, but that can cause issues for
 835                  * r close to 0 for non-interacting pairs.
 836                  */
 837                 real rs, frac, f_lr;
 838                 int  ri;
 839
 840                 rs     = rsq*rinv*ewtabscale;
 841                 ri     = (int)rs;
 842                 frac   = rs - ri;
 843                 f_lr   = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1];
 844                 /* TODO: Currently the Ewald LJ table does not contain
 845                  * the factor 1/6, we should add this.
 846                  */
 847                 FF     = f_lr*rinv/six;
 848                 VV     = (tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr))/six;
 849
 850                 if (ii == jnr)
 851                 {
 852                     /* If we get here, the i particle (ii) has itself (jnr)
 853                      * in its neighborlist. This can only happen with the Verlet
 854                      * scheme, and corresponds to a self-interaction that will
 855                      * occur twice. Scale it down by 50% to only include it once.
 856                      */
 857                     VV *= half;
 858                 }
 859
 860                 for (i = 0; i < NSTATES; i++)
 861                 {
 862                     c6grid      = nbfp_grid[tj[i]];
 863                     vvtot      += LFV[i]*c6grid*VV;
 864                     Fscal      += LFV[i]*c6grid*FF;
 865                     dvdl_vdw   += (DLF[i]*c6grid)*VV;
 866                 }
 867             }
 868
 869             if (bDoForces)
 870             {
 871                 tx         = Fscal*dx;
 872                 ty         = Fscal*dy;
 873                 tz         = Fscal*dz;
 874                 fix        = fix + tx;
 875                 fiy        = fiy + ty;
 876                 fiz        = fiz + tz;
 877                 /* OpenMP atomics are expensive, but this kernels is also
 878                  * expensive, so we can take this hit, instead of using
 879                  * thread-local output buffers and extra reduction.
 880                  */
 881 #pragma omp atomic
 882                 f[j3]     -= tx;
 883 #pragma omp atomic
 884                 f[j3+1]   -= ty;
 885 #pragma omp atomic
 886                 f[j3+2]   -= tz;
 887             }
 888         }
 889
 890         /* The atomics below are expensive with many OpenMP threads.
 891          * Here unperturbed i-particles will usually only have a few
 892          * (perturbed) j-particles in the list. Thus with a buffered list
 893          * we can skip a significant number of i-reductions with a check.
 894          */
 895         if (npair_within_cutoff > 0)
 896         {
 897             if (bDoForces)
 898             {
 899 #pragma omp atomic
 900                 f[ii3]        += fix;
 901 #pragma omp atomic
 902                 f[ii3+1]      += fiy;
 903 #pragma omp atomic
 904                 f[ii3+2]      += fiz;
 905             }
 906             if (bDoShiftForces)
 907             {
 908 #pragma omp atomic
 909                 fshift[is3]   += fix;
 910 #pragma omp atomic
 911                 fshift[is3+1] += fiy;
 912 #pragma omp atomic
 913                 fshift[is3+2] += fiz;
 914             }
 915             if (bDoPotential)
 916             {
 917                 ggid               = gid[n];
 918 #pragma omp atomic
 919                 Vc[ggid]          += vctot;
 920 #pragma omp atomic
 921                 Vv[ggid]          += vvtot;
 922             }
 923         }
 924     }
 925
 926 #pragma omp atomic
 927     dvdl[efptCOUL]     += dvdl_coul;
 928  #pragma omp atomic
 929     dvdl[efptVDW]      += dvdl_vdw;
 930
 931     /* Estimate flops, average for free energy stuff:
 932      * 12  flops per outer iteration
 933      * 150 flops per inner iteration
 934      */
 935 #pragma omp atomic
 936     inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150);
 937 }
 938
 939 real
 940 nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, real alpha_vdw,
 941                                real tabscale, real *vftab,
 942                                real qqA, real c6A, real c12A, real qqB, real c6B, real c12B,
 943                                real LFC[2], real LFV[2], real DLF[2],
 944                                real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], real dlfac_vdw[2],
 945                                real sigma6_def, real sigma6_min, real sigma2_def, real sigma2_min,
 946                                real *velectot, real *vvdwtot, real *dvdl)
 947 {
 948     real       r, rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
 949     real       qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2], sigma_powm2[2];
 950     real       alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
 951     real       rpinv, r_coul, r_vdw, velecsum, vvdwsum;
 952     real       fscal_vdw[2], fscal_elec[2];
 953     real       velec[2], vvdw[2];
 954     int        i, ntab;
 955
 956     const real half        = 0.5;
 957     const real one         = 1.0;
 958     const real two         = 2.0;
 959     const real six         = 6.0;
 960     const real fourtyeight = 48.0;
 961
 962     qq[0]    = qqA;
 963     qq[1]    = qqB;
 964     c6[0]    = c6A;
 965     c6[1]    = c6B;
 966     c12[0]   = c12A;
 967     c12[1]   = c12B;
 968
 969     if (sc_r_power == six)
 970     {
 971         rpm2             = r2*r2;   /* r4 */
 972         rp               = rpm2*r2; /* r6 */
 973     }
 974     else if (sc_r_power == fourtyeight)
 975     {
 976         rp               = r2*r2*r2; /* r6 */
 977         rp               = rp*rp;    /* r12 */
 978         rp               = rp*rp;    /* r24 */
 979         rp               = rp*rp;    /* r48 */
 980         rpm2             = rp/r2;    /* r46 */
 981     }
 982     else
 983     {
 984         rp             = pow(r2, half*sc_r_power);  /* not currently supported as input, but can handle it */
 985         rpm2           = rp/r2;
 986     }
 987
 988     /* Loop over state A(0) and B(1) */
 989     for (i = 0; i < 2; i++)
 990     {
 991         if ((c6[i] > 0) && (c12[i] > 0))
 992         {
 993             /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
 994              * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
 995              */
 996             sigma6[i]       = half*c12[i]/c6[i];
 997             sigma2[i]       = pow(half*c12[i]/c6[i], 1.0/3.0);
 998             /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 999                what data to store externally.  Can't be fixed without larger scale changes, so not 5.0 */
1000             if (sigma6[i] < sigma6_min)   /* for disappearing coul and vdw with soft core at the same time */
1001             {
1002                 sigma6[i] = sigma6_min;
1003                 sigma2[i] = sigma2_min;
1004             }
1005         }
1006         else
1007         {
1008             sigma6[i]       = sigma6_def;
1009             sigma2[i]       = sigma2_def;
1010         }
1011         if (sc_r_power == six)
1012         {
1013             sigma_pow[i]    = sigma6[i];
1014             sigma_powm2[i]  = sigma6[i]/sigma2[i];
1015         }
1016         else if (sc_r_power == fourtyeight)
1017         {
1018             sigma_pow[i]    = sigma6[i]*sigma6[i];       /* sigma^12 */
1019             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
1020             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
1021             sigma_powm2[i]  = sigma_pow[i]/sigma2[i];
1022         }
1023         else
1024         {    /* not really supported as input, but in here for testing the general case*/
1025             sigma_pow[i]    = pow(sigma2[i], sc_r_power/2);
1026             sigma_powm2[i]  = sigma_pow[i]/(sigma2[i]);
1027         }
1028     }
1029
1030     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
1031     if ((c12[0] > 0) && (c12[1] > 0))
1032     {
1033         alpha_vdw_eff    = 0;
1034         alpha_coul_eff   = 0;
1035     }
1036     else
1037     {
1038         alpha_vdw_eff    = alpha_vdw;
1039         alpha_coul_eff   = alpha_coul;
1040     }
1041
1042     /* Loop over A and B states again */
1043     for (i = 0; i < 2; i++)
1044     {
1045         fscal_elec[i] = 0;
1046         fscal_vdw[i]  = 0;
1047         velec[i]      = 0;
1048         vvdw[i]       = 0;
1049
1050         /* Only spend time on A or B state if it is non-zero */
1051         if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
1052         {
1053             /* Coulomb */
1054             rpinv            = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
1055             r_coul           = pow(rpinv, -one/sc_r_power);
1056
1057             /* Electrostatics table lookup data */
1058             rtab             = r_coul*tabscale;
1059             ntab             = rtab;
1060             eps              = rtab-ntab;
1061             eps2             = eps*eps;
1062             ntab             = 12*ntab;
1063             /* Electrostatics */
1064             Y                = vftab[ntab];
1065             F                = vftab[ntab+1];
1066             Geps             = eps*vftab[ntab+2];
1067             Heps2            = eps2*vftab[ntab+3];
1068             Fp               = F+Geps+Heps2;
1069             VV               = Y+eps*Fp;
1070             FF               = Fp+Geps+two*Heps2;
1071             velec[i]         = qq[i]*VV;
1072             fscal_elec[i]    = -qq[i]*FF*r_coul*rpinv*tabscale;
1073
1074             /* Vdw */
1075             rpinv            = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
1076             r_vdw            = pow(rpinv, -one/sc_r_power);
1077             /* Vdw table lookup data */
1078             rtab             = r_vdw*tabscale;
1079             ntab             = rtab;
1080             eps              = rtab-ntab;
1081             eps2             = eps*eps;
1082             ntab             = 12*ntab;
1083             /* Dispersion */
1084             Y                = vftab[ntab+4];
1085             F                = vftab[ntab+5];
1086             Geps             = eps*vftab[ntab+6];
1087             Heps2            = eps2*vftab[ntab+7];
1088             Fp               = F+Geps+Heps2;
1089             VV               = Y+eps*Fp;
1090             FF               = Fp+Geps+two*Heps2;
1091             vvdw[i]          = c6[i]*VV;
1092             fscal_vdw[i]     = -c6[i]*FF;
1093
1094             /* Repulsion */
1095             Y                = vftab[ntab+8];
1096             F                = vftab[ntab+9];
1097             Geps             = eps*vftab[ntab+10];
1098             Heps2            = eps2*vftab[ntab+11];
1099             Fp               = F+Geps+Heps2;
1100             VV               = Y+eps*Fp;
1101             FF               = Fp+Geps+two*Heps2;
1102             vvdw[i]         += c12[i]*VV;
1103             fscal_vdw[i]    -= c12[i]*FF;
1104             fscal_vdw[i]    *= r_vdw*rpinv*tabscale;
1105         }
1106     }
1107     /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
1108     /* Assemble A and B states */
1109     velecsum  = 0;
1110     vvdwsum   = 0;
1111     dvdl_coul = 0;
1112     dvdl_vdw  = 0;
1113     fscal     = 0;
1114     for (i = 0; i < 2; i++)
1115     {
1116         velecsum      += LFC[i]*velec[i];
1117         vvdwsum       += LFV[i]*vvdw[i];
1118
1119         fscal         += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2;
1120
1121         dvdl_coul     += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i];
1122         dvdl_vdw      += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i];
1123     }
1124
1125     dvdl[efptCOUL]     += dvdl_coul;
1126     dvdl[efptVDW]      += dvdl_vdw;
1127
1128     *velectot           = velecsum;
1129     *vvdwtot            = vvdwsum;
1130
1131     return fscal;
1132 }