src/gromacs/mdlib/forcerec.cpp

   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,2015, 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 "config.h"
  40
  41 #include <assert.h>
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "gromacs/domdec/domdec.h"
  46 #include "gromacs/ewald/ewald.h"
  47 #include "gromacs/gmxlib/cuda_tools/pmalloc_cuda.h"
  48 #include "gromacs/gmxlib/gpu_utils/gpu_utils.h"
  49 #include "gromacs/legacyheaders/copyrite.h"
  50 #include "gromacs/legacyheaders/force.h"
  51 #include "gromacs/legacyheaders/gmx_detect_hardware.h"
  52 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  53 #include "gromacs/legacyheaders/inputrec.h"
  54 #include "gromacs/legacyheaders/macros.h"
  55 #include "gromacs/legacyheaders/md_logging.h"
  56 #include "gromacs/legacyheaders/md_support.h"
  57 #include "gromacs/legacyheaders/names.h"
  58 #include "gromacs/legacyheaders/network.h"
  59 #include "gromacs/legacyheaders/nonbonded.h"
  60 #include "gromacs/legacyheaders/ns.h"
  61 #include "gromacs/legacyheaders/qmmm.h"
  62 #include "gromacs/legacyheaders/tables.h"
  63 #include "gromacs/legacyheaders/txtdump.h"
  64 #include "gromacs/legacyheaders/typedefs.h"
  65 #include "gromacs/legacyheaders/types/commrec.h"
  66 #include "gromacs/legacyheaders/types/nbnxn_cuda_types_ext.h"
  67 #include "gromacs/math/calculate-ewald-splitting-coefficient.h"
  68 #include "gromacs/math/units.h"
  69 #include "gromacs/math/utilities.h"
  70 #include "gromacs/math/vec.h"
  71 #include "gromacs/mdlib/forcerec-threading.h"
  72 #include "gromacs/mdlib/nb_verlet.h"
  73 #include "gromacs/mdlib/nbnxn_atomdata.h"
  74 #include "gromacs/mdlib/nbnxn_consts.h"
  75 #include "gromacs/mdlib/nbnxn_search.h"
  76 #include "gromacs/mdlib/nbnxn_simd.h"
  77 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.h"
  78 #include "gromacs/pbcutil/ishift.h"
  79 #include "gromacs/pbcutil/pbc.h"
  80 #include "gromacs/simd/simd.h"
  81 #include "gromacs/topology/mtop_util.h"
  82 #include "gromacs/utility/fatalerror.h"
  83 #include "gromacs/utility/smalloc.h"
  84
  85 t_forcerec *mk_forcerec(void)
  86 {
  87     t_forcerec *fr;
  88
  89     snew(fr, 1);
  90
  91     return fr;
  92 }
  93
  94 #ifdef DEBUG
  95 static void pr_nbfp(FILE *fp, real *nbfp, gmx_bool bBHAM, int atnr)
  96 {
  97     int i, j;
  98
  99     for (i = 0; (i < atnr); i++)
 100     {
 101         for (j = 0; (j < atnr); j++)
 102         {
 103             fprintf(fp, "%2d - %2d", i, j);
 104             if (bBHAM)
 105             {
 106                 fprintf(fp, "  a=%10g, b=%10g, c=%10g\n", BHAMA(nbfp, atnr, i, j),
 107                         BHAMB(nbfp, atnr, i, j), BHAMC(nbfp, atnr, i, j)/6.0);
 108             }
 109             else
 110             {
 111                 fprintf(fp, "  c6=%10g, c12=%10g\n", C6(nbfp, atnr, i, j)/6.0,
 112                         C12(nbfp, atnr, i, j)/12.0);
 113             }
 114         }
 115     }
 116 }
 117 #endif
 118
 119 static real *mk_nbfp(const gmx_ffparams_t *idef, gmx_bool bBHAM)
 120 {
 121     real *nbfp;
 122     int   i, j, k, atnr;
 123
 124     atnr = idef->atnr;
 125     if (bBHAM)
 126     {
 127         snew(nbfp, 3*atnr*atnr);
 128         for (i = k = 0; (i < atnr); i++)
 129         {
 130             for (j = 0; (j < atnr); j++, k++)
 131             {
 132                 BHAMA(nbfp, atnr, i, j) = idef->iparams[k].bham.a;
 133                 BHAMB(nbfp, atnr, i, j) = idef->iparams[k].bham.b;
 134                 /* nbfp now includes the 6.0 derivative prefactor */
 135                 BHAMC(nbfp, atnr, i, j) = idef->iparams[k].bham.c*6.0;
 136             }
 137         }
 138     }
 139     else
 140     {
 141         snew(nbfp, 2*atnr*atnr);
 142         for (i = k = 0; (i < atnr); i++)
 143         {
 144             for (j = 0; (j < atnr); j++, k++)
 145             {
 146                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
 147                 C6(nbfp, atnr, i, j)   = idef->iparams[k].lj.c6*6.0;
 148                 C12(nbfp, atnr, i, j)  = idef->iparams[k].lj.c12*12.0;
 149             }
 150         }
 151     }
 152
 153     return nbfp;
 154 }
 155
 156 static real *make_ljpme_c6grid(const gmx_ffparams_t *idef, t_forcerec *fr)
 157 {
 158     int        i, j, k, atnr;
 159     real       c6, c6i, c6j, c12i, c12j, epsi, epsj, sigmai, sigmaj;
 160     real      *grid;
 161     const real oneOverSix = 1.0 / 6.0;
 162
 163     /* For LJ-PME simulations, we correct the energies with the reciprocal space
 164      * inside of the cut-off. To do this the non-bonded kernels needs to have
 165      * access to the C6-values used on the reciprocal grid in pme.c
 166      */
 167
 168     atnr = idef->atnr;
 169     snew(grid, 2*atnr*atnr);
 170     for (i = k = 0; (i < atnr); i++)
 171     {
 172         for (j = 0; (j < atnr); j++, k++)
 173         {
 174             c6i  = idef->iparams[i*(atnr+1)].lj.c6;
 175             c12i = idef->iparams[i*(atnr+1)].lj.c12;
 176             c6j  = idef->iparams[j*(atnr+1)].lj.c6;
 177             c12j = idef->iparams[j*(atnr+1)].lj.c12;
 178             c6   = sqrt(c6i * c6j);
 179             if (fr->ljpme_combination_rule == eljpmeLB
 180                 && !gmx_numzero(c6) && !gmx_numzero(c12i) && !gmx_numzero(c12j))
 181             {
 182                 sigmai = pow(c12i / c6i, oneOverSix);
 183                 sigmaj = pow(c12j / c6j, oneOverSix);
 184                 epsi   = c6i * c6i / c12i;
 185                 epsj   = c6j * c6j / c12j;
 186                 c6     = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 6);
 187             }
 188             /* Store the elements at the same relative positions as C6 in nbfp in order
 189              * to simplify access in the kernels
 190              */
 191             grid[2*(atnr*i+j)] = c6*6.0;
 192         }
 193     }
 194     return grid;
 195 }
 196
 197 static real *mk_nbfp_combination_rule(const gmx_ffparams_t *idef, int comb_rule)
 198 {
 199     real      *nbfp;
 200     int        i, j, atnr;
 201     real       c6i, c6j, c12i, c12j, epsi, epsj, sigmai, sigmaj;
 202     real       c6, c12;
 203     const real oneOverSix = 1.0 / 6.0;
 204
 205     atnr = idef->atnr;
 206     snew(nbfp, 2*atnr*atnr);
 207     for (i = 0; i < atnr; ++i)
 208     {
 209         for (j = 0; j < atnr; ++j)
 210         {
 211             c6i  = idef->iparams[i*(atnr+1)].lj.c6;
 212             c12i = idef->iparams[i*(atnr+1)].lj.c12;
 213             c6j  = idef->iparams[j*(atnr+1)].lj.c6;
 214             c12j = idef->iparams[j*(atnr+1)].lj.c12;
 215             c6   = sqrt(c6i  * c6j);
 216             c12  = sqrt(c12i * c12j);
 217             if (comb_rule == eCOMB_ARITHMETIC
 218                 && !gmx_numzero(c6) && !gmx_numzero(c12))
 219             {
 220                 sigmai = pow(c12i / c6i, oneOverSix);
 221                 sigmaj = pow(c12j / c6j, oneOverSix);
 222                 epsi   = c6i * c6i / c12i;
 223                 epsj   = c6j * c6j / c12j;
 224                 c6     = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 6);
 225                 c12    = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 12);
 226             }
 227             C6(nbfp, atnr, i, j)   = c6*6.0;
 228             C12(nbfp, atnr, i, j)  = c12*12.0;
 229         }
 230     }
 231     return nbfp;
 232 }
 233
 234 /* This routine sets fr->solvent_opt to the most common solvent in the
 235  * system, e.g. esolSPC or esolTIP4P. It will also mark each charge group in
 236  * the fr->solvent_type array with the correct type (or esolNO).
 237  *
 238  * Charge groups that fulfill the conditions but are not identical to the
 239  * most common one will be marked as esolNO in the solvent_type array.
 240  *
 241  * TIP3p is identical to SPC for these purposes, so we call it
 242  * SPC in the arrays (Apologies to Bill Jorgensen ;-)
 243  *
 244  * NOTE: QM particle should not
 245  * become an optimized solvent. Not even if there is only one charge
 246  * group in the Qm
 247  */
 248
 249 typedef struct
 250 {
 251     int    model;
 252     int    count;
 253     int    vdwtype[4];
 254     real   charge[4];
 255 } solvent_parameters_t;
 256
 257 static void
 258 check_solvent_cg(const gmx_moltype_t    *molt,
 259                  int                     cg0,
 260                  int                     nmol,
 261                  const unsigned char    *qm_grpnr,
 262                  const t_grps           *qm_grps,
 263                  t_forcerec   *          fr,
 264                  int                    *n_solvent_parameters,
 265                  solvent_parameters_t  **solvent_parameters_p,
 266                  int                     cginfo,
 267                  int                    *cg_sp)
 268 {
 269     t_atom               *atom;
 270     int                   j, k;
 271     int                   j0, j1, nj;
 272     gmx_bool              perturbed;
 273     gmx_bool              has_vdw[4];
 274     gmx_bool              match;
 275     real                  tmp_charge[4]  = { 0.0 }; /* init to zero to make gcc4.8 happy */
 276     int                   tmp_vdwtype[4] = { 0 };   /* init to zero to make gcc4.8 happy */
 277     int                   tjA;
 278     gmx_bool              qm;
 279     solvent_parameters_t *solvent_parameters;
 280
 281     /* We use a list with parameters for each solvent type.
 282      * Every time we discover a new molecule that fulfills the basic
 283      * conditions for a solvent we compare with the previous entries
 284      * in these lists. If the parameters are the same we just increment
 285      * the counter for that type, and otherwise we create a new type
 286      * based on the current molecule.
 287      *
 288      * Once we've finished going through all molecules we check which
 289      * solvent is most common, and mark all those molecules while we
 290      * clear the flag on all others.
 291      */
 292
 293     solvent_parameters = *solvent_parameters_p;
 294
 295     /* Mark the cg first as non optimized */
 296     *cg_sp = -1;
 297
 298     /* Check if this cg has no exclusions with atoms in other charge groups
 299      * and all atoms inside the charge group excluded.
 300      * We only have 3 or 4 atom solvent loops.
 301      */
 302     if (GET_CGINFO_EXCL_INTER(cginfo) ||
 303         !GET_CGINFO_EXCL_INTRA(cginfo))
 304     {
 305         return;
 306     }
 307
 308     /* Get the indices of the first atom in this charge group */
 309     j0     = molt->cgs.index[cg0];
 310     j1     = molt->cgs.index[cg0+1];
 311
 312     /* Number of atoms in our molecule */
 313     nj     = j1 - j0;
 314
 315     if (debug)
 316     {
 317         fprintf(debug,
 318                 "Moltype '%s': there are %d atoms in this charge group\n",
 319                 *molt->name, nj);
 320     }
 321
 322     /* Check if it could be an SPC (3 atoms) or TIP4p (4) water,
 323      * otherwise skip it.
 324      */
 325     if (nj < 3 || nj > 4)
 326     {
 327         return;
 328     }
 329
 330     /* Check if we are doing QM on this group */
 331     qm = FALSE;
 332     if (qm_grpnr != NULL)
 333     {
 334         for (j = j0; j < j1 && !qm; j++)
 335         {
 336             qm = (qm_grpnr[j] < qm_grps->nr - 1);
 337         }
 338     }
 339     /* Cannot use solvent optimization with QM */
 340     if (qm)
 341     {
 342         return;
 343     }
 344
 345     atom = molt->atoms.atom;
 346
 347     /* Still looks like a solvent, time to check parameters */
 348
 349     /* If it is perturbed (free energy) we can't use the solvent loops,
 350      * so then we just skip to the next molecule.
 351      */
 352     perturbed = FALSE;
 353
 354     for (j = j0; j < j1 && !perturbed; j++)
 355     {
 356         perturbed = PERTURBED(atom[j]);
 357     }
 358
 359     if (perturbed)
 360     {
 361         return;
 362     }
 363
 364     /* Now it's only a question if the VdW and charge parameters
 365      * are OK. Before doing the check we compare and see if they are
 366      * identical to a possible previous solvent type.
 367      * First we assign the current types and charges.
 368      */
 369     for (j = 0; j < nj; j++)
 370     {
 371         tmp_vdwtype[j] = atom[j0+j].type;
 372         tmp_charge[j]  = atom[j0+j].q;
 373     }
 374
 375     /* Does it match any previous solvent type? */
 376     for (k = 0; k < *n_solvent_parameters; k++)
 377     {
 378         match = TRUE;
 379
 380
 381         /* We can only match SPC with 3 atoms and TIP4p with 4 atoms */
 382         if ( (solvent_parameters[k].model == esolSPC   && nj != 3)  ||
 383              (solvent_parameters[k].model == esolTIP4P && nj != 4) )
 384         {
 385             match = FALSE;
 386         }
 387
 388         /* Check that types & charges match for all atoms in molecule */
 389         for (j = 0; j < nj && match == TRUE; j++)
 390         {
 391             if (tmp_vdwtype[j] != solvent_parameters[k].vdwtype[j])
 392             {
 393                 match = FALSE;
 394             }
 395             if (tmp_charge[j] != solvent_parameters[k].charge[j])
 396             {
 397                 match = FALSE;
 398             }
 399         }
 400         if (match == TRUE)
 401         {
 402             /* Congratulations! We have a matched solvent.
 403              * Flag it with this type for later processing.
 404              */
 405             *cg_sp = k;
 406             solvent_parameters[k].count += nmol;
 407
 408             /* We are done with this charge group */
 409             return;
 410         }
 411     }
 412
 413     /* If we get here, we have a tentative new solvent type.
 414      * Before we add it we must check that it fulfills the requirements
 415      * of the solvent optimized loops. First determine which atoms have
 416      * VdW interactions.
 417      */
 418     for (j = 0; j < nj; j++)
 419     {
 420         has_vdw[j] = FALSE;
 421         tjA        = tmp_vdwtype[j];
 422
 423         /* Go through all other tpes and see if any have non-zero
 424          * VdW parameters when combined with this one.
 425          */
 426         for (k = 0; k < fr->ntype && (has_vdw[j] == FALSE); k++)
 427         {
 428             /* We already checked that the atoms weren't perturbed,
 429              * so we only need to check state A now.
 430              */
 431             if (fr->bBHAM)
 432             {
 433                 has_vdw[j] = (has_vdw[j] ||
 434                               (BHAMA(fr->nbfp, fr->ntype, tjA, k) != 0.0) ||
 435                               (BHAMB(fr->nbfp, fr->ntype, tjA, k) != 0.0) ||
 436                               (BHAMC(fr->nbfp, fr->ntype, tjA, k) != 0.0));
 437             }
 438             else
 439             {
 440                 /* Standard LJ */
 441                 has_vdw[j] = (has_vdw[j] ||
 442                               (C6(fr->nbfp, fr->ntype, tjA, k)  != 0.0) ||
 443                               (C12(fr->nbfp, fr->ntype, tjA, k) != 0.0));
 444             }
 445         }
 446     }
 447
 448     /* Now we know all we need to make the final check and assignment. */
 449     if (nj == 3)
 450     {
 451         /* So, is it an SPC?
 452          * For this we require thatn all atoms have charge,
 453          * the charges on atom 2 & 3 should be the same, and only
 454          * atom 1 might have VdW.
 455          */
 456         if (has_vdw[1] == FALSE &&
 457             has_vdw[2] == FALSE &&
 458             tmp_charge[0]  != 0 &&
 459             tmp_charge[1]  != 0 &&
 460             tmp_charge[2]  == tmp_charge[1])
 461         {
 462             srenew(solvent_parameters, *n_solvent_parameters+1);
 463             solvent_parameters[*n_solvent_parameters].model = esolSPC;
 464             solvent_parameters[*n_solvent_parameters].count = nmol;
 465             for (k = 0; k < 3; k++)
 466             {
 467                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 468                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 469             }
 470
 471             *cg_sp = *n_solvent_parameters;
 472             (*n_solvent_parameters)++;
 473         }
 474     }
 475     else if (nj == 4)
 476     {
 477         /* Or could it be a TIP4P?
 478          * For this we require thatn atoms 2,3,4 have charge, but not atom 1.
 479          * Only atom 1 mght have VdW.
 480          */
 481         if (has_vdw[1] == FALSE &&
 482             has_vdw[2] == FALSE &&
 483             has_vdw[3] == FALSE &&
 484             tmp_charge[0]  == 0 &&
 485             tmp_charge[1]  != 0 &&
 486             tmp_charge[2]  == tmp_charge[1] &&
 487             tmp_charge[3]  != 0)
 488         {
 489             srenew(solvent_parameters, *n_solvent_parameters+1);
 490             solvent_parameters[*n_solvent_parameters].model = esolTIP4P;
 491             solvent_parameters[*n_solvent_parameters].count = nmol;
 492             for (k = 0; k < 4; k++)
 493             {
 494                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 495                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 496             }
 497
 498             *cg_sp = *n_solvent_parameters;
 499             (*n_solvent_parameters)++;
 500         }
 501     }
 502
 503     *solvent_parameters_p = solvent_parameters;
 504 }
 505
 506 static void
 507 check_solvent(FILE  *                fp,
 508               const gmx_mtop_t  *    mtop,
 509               t_forcerec  *          fr,
 510               cginfo_mb_t           *cginfo_mb)
 511 {
 512     const t_block     *   cgs;
 513     const gmx_moltype_t  *molt;
 514     int                   mb, mol, cg_mol, at_offset, am, cgm, i, nmol_ch, nmol;
 515     int                   n_solvent_parameters;
 516     solvent_parameters_t *solvent_parameters;
 517     int                 **cg_sp;
 518     int                   bestsp, bestsol;
 519
 520     if (debug)
 521     {
 522         fprintf(debug, "Going to determine what solvent types we have.\n");
 523     }
 524
 525     n_solvent_parameters = 0;
 526     solvent_parameters   = NULL;
 527     /* Allocate temporary array for solvent type */
 528     snew(cg_sp, mtop->nmolblock);
 529
 530     at_offset = 0;
 531     for (mb = 0; mb < mtop->nmolblock; mb++)
 532     {
 533         molt = &mtop->moltype[mtop->molblock[mb].type];
 534         cgs  = &molt->cgs;
 535         /* Here we have to loop over all individual molecules
 536          * because we need to check for QMMM particles.
 537          */
 538         snew(cg_sp[mb], cginfo_mb[mb].cg_mod);
 539         nmol_ch = cginfo_mb[mb].cg_mod/cgs->nr;
 540         nmol    = mtop->molblock[mb].nmol/nmol_ch;
 541         for (mol = 0; mol < nmol_ch; mol++)
 542         {
 543             cgm = mol*cgs->nr;
 544             am  = mol*cgs->index[cgs->nr];
 545             for (cg_mol = 0; cg_mol < cgs->nr; cg_mol++)
 546             {
 547                 check_solvent_cg(molt, cg_mol, nmol,
 548                                  mtop->groups.grpnr[egcQMMM] ?
 549                                  mtop->groups.grpnr[egcQMMM]+at_offset+am : 0,
 550                                  &mtop->groups.grps[egcQMMM],
 551                                  fr,
 552                                  &n_solvent_parameters, &solvent_parameters,
 553                                  cginfo_mb[mb].cginfo[cgm+cg_mol],
 554                                  &cg_sp[mb][cgm+cg_mol]);
 555             }
 556         }
 557         at_offset += cgs->index[cgs->nr];
 558     }
 559
 560     /* Puh! We finished going through all charge groups.
 561      * Now find the most common solvent model.
 562      */
 563
 564     /* Most common solvent this far */
 565     bestsp = -2;
 566     for (i = 0; i < n_solvent_parameters; i++)
 567     {
 568         if (bestsp == -2 ||
 569             solvent_parameters[i].count > solvent_parameters[bestsp].count)
 570         {
 571             bestsp = i;
 572         }
 573     }
 574
 575     if (bestsp >= 0)
 576     {
 577         bestsol = solvent_parameters[bestsp].model;
 578     }
 579     else
 580     {
 581         bestsol = esolNO;
 582     }
 583
 584     fr->nWatMol = 0;
 585     for (mb = 0; mb < mtop->nmolblock; mb++)
 586     {
 587         cgs  = &mtop->moltype[mtop->molblock[mb].type].cgs;
 588         nmol = (mtop->molblock[mb].nmol*cgs->nr)/cginfo_mb[mb].cg_mod;
 589         for (i = 0; i < cginfo_mb[mb].cg_mod; i++)
 590         {
 591             if (cg_sp[mb][i] == bestsp)
 592             {
 593                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i], bestsol);
 594                 fr->nWatMol += nmol;
 595             }
 596             else
 597             {
 598                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i], esolNO);
 599             }
 600         }
 601         sfree(cg_sp[mb]);
 602     }
 603     sfree(cg_sp);
 604
 605     if (bestsol != esolNO && fp != NULL)
 606     {
 607         fprintf(fp, "\nEnabling %s-like water optimization for %d molecules.\n\n",
 608                 esol_names[bestsol],
 609                 solvent_parameters[bestsp].count);
 610     }
 611
 612     sfree(solvent_parameters);
 613     fr->solvent_opt = bestsol;
 614 }
 615
 616 enum {
 617     acNONE = 0, acCONSTRAINT, acSETTLE
 618 };
 619
 620 static cginfo_mb_t *init_cginfo_mb(FILE *fplog, const gmx_mtop_t *mtop,
 621                                    t_forcerec *fr, gmx_bool bNoSolvOpt,
 622                                    gmx_bool *bFEP_NonBonded,
 623                                    gmx_bool *bExcl_IntraCGAll_InterCGNone)
 624 {
 625     const t_block        *cgs;
 626     const t_blocka       *excl;
 627     const gmx_moltype_t  *molt;
 628     const gmx_molblock_t *molb;
 629     cginfo_mb_t          *cginfo_mb;
 630     gmx_bool             *type_VDW;
 631     int                  *cginfo;
 632     int                   cg_offset, a_offset;
 633     int                   mb, m, cg, a0, a1, gid, ai, j, aj, excl_nalloc;
 634     int                  *a_con;
 635     int                   ftype;
 636     int                   ia;
 637     gmx_bool              bId, *bExcl, bExclIntraAll, bExclInter, bHaveVDW, bHaveQ, bHavePerturbedAtoms;
 638
 639     snew(cginfo_mb, mtop->nmolblock);
 640
 641     snew(type_VDW, fr->ntype);
 642     for (ai = 0; ai < fr->ntype; ai++)
 643     {
 644         type_VDW[ai] = FALSE;
 645         for (j = 0; j < fr->ntype; j++)
 646         {
 647             type_VDW[ai] = type_VDW[ai] ||
 648                 fr->bBHAM ||
 649                 C6(fr->nbfp, fr->ntype, ai, j) != 0 ||
 650                 C12(fr->nbfp, fr->ntype, ai, j) != 0;
 651         }
 652     }
 653
 654     *bFEP_NonBonded               = FALSE;
 655     *bExcl_IntraCGAll_InterCGNone = TRUE;
 656
 657     excl_nalloc = 10;
 658     snew(bExcl, excl_nalloc);
 659     cg_offset = 0;
 660     a_offset  = 0;
 661     for (mb = 0; mb < mtop->nmolblock; mb++)
 662     {
 663         molb = &mtop->molblock[mb];
 664         molt = &mtop->moltype[molb->type];
 665         cgs  = &molt->cgs;
 666         excl = &molt->excls;
 667
 668         /* Check if the cginfo is identical for all molecules in this block.
 669          * If so, we only need an array of the size of one molecule.
 670          * Otherwise we make an array of #mol times #cgs per molecule.
 671          */
 672         bId = TRUE;
 673         for (m = 0; m < molb->nmol; m++)
 674         {
 675             int am = m*cgs->index[cgs->nr];
 676             for (cg = 0; cg < cgs->nr; cg++)
 677             {
 678                 a0 = cgs->index[cg];
 679                 a1 = cgs->index[cg+1];
 680                 if (ggrpnr(&mtop->groups, egcENER, a_offset+am+a0) !=
 681                     ggrpnr(&mtop->groups, egcENER, a_offset   +a0))
 682                 {
 683                     bId = FALSE;
 684                 }
 685                 if (mtop->groups.grpnr[egcQMMM] != NULL)
 686                 {
 687                     for (ai = a0; ai < a1; ai++)
 688                     {
 689                         if (mtop->groups.grpnr[egcQMMM][a_offset+am+ai] !=
 690                             mtop->groups.grpnr[egcQMMM][a_offset   +ai])
 691                         {
 692                             bId = FALSE;
 693                         }
 694                     }
 695                 }
 696             }
 697         }
 698
 699         cginfo_mb[mb].cg_start = cg_offset;
 700         cginfo_mb[mb].cg_end   = cg_offset + molb->nmol*cgs->nr;
 701         cginfo_mb[mb].cg_mod   = (bId ? 1 : molb->nmol)*cgs->nr;
 702         snew(cginfo_mb[mb].cginfo, cginfo_mb[mb].cg_mod);
 703         cginfo = cginfo_mb[mb].cginfo;
 704
 705         /* Set constraints flags for constrained atoms */
 706         snew(a_con, molt->atoms.nr);
 707         for (ftype = 0; ftype < F_NRE; ftype++)
 708         {
 709             if (interaction_function[ftype].flags & IF_CONSTRAINT)
 710             {
 711                 int nral;
 712
 713                 nral = NRAL(ftype);
 714                 for (ia = 0; ia < molt->ilist[ftype].nr; ia += 1+nral)
 715                 {
 716                     int a;
 717
 718                     for (a = 0; a < nral; a++)
 719                     {
 720                         a_con[molt->ilist[ftype].iatoms[ia+1+a]] =
 721                             (ftype == F_SETTLE ? acSETTLE : acCONSTRAINT);
 722                     }
 723                 }
 724             }
 725         }
 726
 727         for (m = 0; m < (bId ? 1 : molb->nmol); m++)
 728         {
 729             int cgm = m*cgs->nr;
 730             int am  = m*cgs->index[cgs->nr];
 731             for (cg = 0; cg < cgs->nr; cg++)
 732             {
 733                 a0 = cgs->index[cg];
 734                 a1 = cgs->index[cg+1];
 735
 736                 /* Store the energy group in cginfo */
 737                 gid = ggrpnr(&mtop->groups, egcENER, a_offset+am+a0);
 738                 SET_CGINFO_GID(cginfo[cgm+cg], gid);
 739
 740                 /* Check the intra/inter charge group exclusions */
 741                 if (a1-a0 > excl_nalloc)
 742                 {
 743                     excl_nalloc = a1 - a0;
 744                     srenew(bExcl, excl_nalloc);
 745                 }
 746                 /* bExclIntraAll: all intra cg interactions excluded
 747                  * bExclInter:    any inter cg interactions excluded
 748                  */
 749                 bExclIntraAll       = TRUE;
 750                 bExclInter          = FALSE;
 751                 bHaveVDW            = FALSE;
 752                 bHaveQ              = FALSE;
 753                 bHavePerturbedAtoms = FALSE;
 754                 for (ai = a0; ai < a1; ai++)
 755                 {
 756                     /* Check VDW and electrostatic interactions */
 757                     bHaveVDW = bHaveVDW || (type_VDW[molt->atoms.atom[ai].type] ||
 758                                             type_VDW[molt->atoms.atom[ai].typeB]);
 759                     bHaveQ  = bHaveQ    || (molt->atoms.atom[ai].q != 0 ||
 760                                             molt->atoms.atom[ai].qB != 0);
 761
 762                     bHavePerturbedAtoms = bHavePerturbedAtoms || (PERTURBED(molt->atoms.atom[ai]) != 0);
 763
 764                     /* Clear the exclusion list for atom ai */
 765                     for (aj = a0; aj < a1; aj++)
 766                     {
 767                         bExcl[aj-a0] = FALSE;
 768                     }
 769                     /* Loop over all the exclusions of atom ai */
 770                     for (j = excl->index[ai]; j < excl->index[ai+1]; j++)
 771                     {
 772                         aj = excl->a[j];
 773                         if (aj < a0 || aj >= a1)
 774                         {
 775                             bExclInter = TRUE;
 776                         }
 777                         else
 778                         {
 779                             bExcl[aj-a0] = TRUE;
 780                         }
 781                     }
 782                     /* Check if ai excludes a0 to a1 */
 783                     for (aj = a0; aj < a1; aj++)
 784                     {
 785                         if (!bExcl[aj-a0])
 786                         {
 787                             bExclIntraAll = FALSE;
 788                         }
 789                     }
 790
 791                     switch (a_con[ai])
 792                     {
 793                         case acCONSTRAINT:
 794                             SET_CGINFO_CONSTR(cginfo[cgm+cg]);
 795                             break;
 796                         case acSETTLE:
 797                             SET_CGINFO_SETTLE(cginfo[cgm+cg]);
 798                             break;
 799                         default:
 800                             break;
 801                     }
 802                 }
 803                 if (bExclIntraAll)
 804                 {
 805                     SET_CGINFO_EXCL_INTRA(cginfo[cgm+cg]);
 806                 }
 807                 if (bExclInter)
 808                 {
 809                     SET_CGINFO_EXCL_INTER(cginfo[cgm+cg]);
 810                 }
 811                 if (a1 - a0 > MAX_CHARGEGROUP_SIZE)
 812                 {
 813                     /* The size in cginfo is currently only read with DD */
 814                     gmx_fatal(FARGS, "A charge group has size %d which is larger than the limit of %d atoms", a1-a0, MAX_CHARGEGROUP_SIZE);
 815                 }
 816                 if (bHaveVDW)
 817                 {
 818                     SET_CGINFO_HAS_VDW(cginfo[cgm+cg]);
 819                 }
 820                 if (bHaveQ)
 821                 {
 822                     SET_CGINFO_HAS_Q(cginfo[cgm+cg]);
 823                 }
 824                 if (bHavePerturbedAtoms && fr->efep != efepNO)
 825                 {
 826                     SET_CGINFO_FEP(cginfo[cgm+cg]);
 827                     *bFEP_NonBonded = TRUE;
 828                 }
 829                 /* Store the charge group size */
 830                 SET_CGINFO_NATOMS(cginfo[cgm+cg], a1-a0);
 831
 832                 if (!bExclIntraAll || bExclInter)
 833                 {
 834                     *bExcl_IntraCGAll_InterCGNone = FALSE;
 835                 }
 836             }
 837         }
 838
 839         sfree(a_con);
 840
 841         cg_offset += molb->nmol*cgs->nr;
 842         a_offset  += molb->nmol*cgs->index[cgs->nr];
 843     }
 844     sfree(bExcl);
 845
 846     /* the solvent optimizer is called after the QM is initialized,
 847      * because we don't want to have the QM subsystemto become an
 848      * optimized solvent
 849      */
 850
 851     check_solvent(fplog, mtop, fr, cginfo_mb);
 852
 853     if (getenv("GMX_NO_SOLV_OPT"))
 854     {
 855         if (fplog)
 856         {
 857             fprintf(fplog, "Found environment variable GMX_NO_SOLV_OPT.\n"
 858                     "Disabling all solvent optimization\n");
 859         }
 860         fr->solvent_opt = esolNO;
 861     }
 862     if (bNoSolvOpt)
 863     {
 864         fr->solvent_opt = esolNO;
 865     }
 866     if (!fr->solvent_opt)
 867     {
 868         for (mb = 0; mb < mtop->nmolblock; mb++)
 869         {
 870             for (cg = 0; cg < cginfo_mb[mb].cg_mod; cg++)
 871             {
 872                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[cg], esolNO);
 873             }
 874         }
 875     }
 876
 877     return cginfo_mb;
 878 }
 879
 880 static int *cginfo_expand(int nmb, cginfo_mb_t *cgi_mb)
 881 {
 882     int  ncg, mb, cg;
 883     int *cginfo;
 884
 885     ncg = cgi_mb[nmb-1].cg_end;
 886     snew(cginfo, ncg);
 887     mb = 0;
 888     for (cg = 0; cg < ncg; cg++)
 889     {
 890         while (cg >= cgi_mb[mb].cg_end)
 891         {
 892             mb++;
 893         }
 894         cginfo[cg] =
 895             cgi_mb[mb].cginfo[(cg - cgi_mb[mb].cg_start) % cgi_mb[mb].cg_mod];
 896     }
 897
 898     return cginfo;
 899 }
 900
 901 static void set_chargesum(FILE *log, t_forcerec *fr, const gmx_mtop_t *mtop)
 902 {
 903     /*This now calculates sum for q and c6*/
 904     double         qsum, q2sum, q, c6sum, c6;
 905     int            mb, nmol, i;
 906     const t_atoms *atoms;
 907
 908     qsum   = 0;
 909     q2sum  = 0;
 910     c6sum  = 0;
 911     for (mb = 0; mb < mtop->nmolblock; mb++)
 912     {
 913         nmol  = mtop->molblock[mb].nmol;
 914         atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 915         for (i = 0; i < atoms->nr; i++)
 916         {
 917             q       = atoms->atom[i].q;
 918             qsum   += nmol*q;
 919             q2sum  += nmol*q*q;
 920             c6      = mtop->ffparams.iparams[atoms->atom[i].type*(mtop->ffparams.atnr+1)].lj.c6;
 921             c6sum  += nmol*c6;
 922         }
 923     }
 924     fr->qsum[0]   = qsum;
 925     fr->q2sum[0]  = q2sum;
 926     fr->c6sum[0]  = c6sum;
 927
 928     if (fr->efep != efepNO)
 929     {
 930         qsum   = 0;
 931         q2sum  = 0;
 932         c6sum  = 0;
 933         for (mb = 0; mb < mtop->nmolblock; mb++)
 934         {
 935             nmol  = mtop->molblock[mb].nmol;
 936             atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 937             for (i = 0; i < atoms->nr; i++)
 938             {
 939                 q       = atoms->atom[i].qB;
 940                 qsum   += nmol*q;
 941                 q2sum  += nmol*q*q;
 942                 c6      = mtop->ffparams.iparams[atoms->atom[i].typeB*(mtop->ffparams.atnr+1)].lj.c6;
 943                 c6sum  += nmol*c6;
 944             }
 945             fr->qsum[1]   = qsum;
 946             fr->q2sum[1]  = q2sum;
 947             fr->c6sum[1]  = c6sum;
 948         }
 949     }
 950     else
 951     {
 952         fr->qsum[1]   = fr->qsum[0];
 953         fr->q2sum[1]  = fr->q2sum[0];
 954         fr->c6sum[1]  = fr->c6sum[0];
 955     }
 956     if (log)
 957     {
 958         if (fr->efep == efepNO)
 959         {
 960             fprintf(log, "System total charge: %.3f\n", fr->qsum[0]);
 961         }
 962         else
 963         {
 964             fprintf(log, "System total charge, top. A: %.3f top. B: %.3f\n",
 965                     fr->qsum[0], fr->qsum[1]);
 966         }
 967     }
 968 }
 969
 970 void update_forcerec(t_forcerec *fr, matrix box)
 971 {
 972     if (fr->eeltype == eelGRF)
 973     {
 974         calc_rffac(NULL, fr->eeltype, fr->epsilon_r, fr->epsilon_rf,
 975                    fr->rcoulomb, fr->temp, fr->zsquare, box,
 976                    &fr->kappa, &fr->k_rf, &fr->c_rf);
 977     }
 978 }
 979
 980 void set_avcsixtwelve(FILE *fplog, t_forcerec *fr, const gmx_mtop_t *mtop)
 981 {
 982     const t_atoms  *atoms, *atoms_tpi;
 983     const t_blocka *excl;
 984     int             mb, nmol, nmolc, i, j, tpi, tpj, j1, j2, k, nexcl, q;
 985     gmx_int64_t     npair, npair_ij, tmpi, tmpj;
 986     double          csix, ctwelve;
 987     int             ntp, *typecount;
 988     gmx_bool        bBHAM;
 989     real           *nbfp;
 990     real           *nbfp_comb = NULL;
 991
 992     ntp   = fr->ntype;
 993     bBHAM = fr->bBHAM;
 994     nbfp  = fr->nbfp;
 995
 996     /* For LJ-PME, we want to correct for the difference between the
 997      * actual C6 values and the C6 values used by the LJ-PME based on
 998      * combination rules. */
 999
1000     if (EVDW_PME(fr->vdwtype))
1001     {
1002         nbfp_comb = mk_nbfp_combination_rule(&mtop->ffparams,
1003                                              (fr->ljpme_combination_rule == eljpmeLB) ? eCOMB_ARITHMETIC : eCOMB_GEOMETRIC);
1004         for (tpi = 0; tpi < ntp; ++tpi)
1005         {
1006             for (tpj = 0; tpj < ntp; ++tpj)
1007             {
1008                 C6(nbfp_comb, ntp, tpi, tpj) =
1009                     C6(nbfp, ntp, tpi, tpj) - C6(nbfp_comb, ntp, tpi, tpj);
1010                 C12(nbfp_comb, ntp, tpi, tpj) = C12(nbfp, ntp, tpi, tpj);
1011             }
1012         }
1013         nbfp = nbfp_comb;
1014     }
1015     for (q = 0; q < (fr->efep == efepNO ? 1 : 2); q++)
1016     {
1017         csix    = 0;
1018         ctwelve = 0;
1019         npair   = 0;
1020         nexcl   = 0;
1021         if (!fr->n_tpi)
1022         {
1023             /* Count the types so we avoid natoms^2 operations */
1024             snew(typecount, ntp);
1025             gmx_mtop_count_atomtypes(mtop, q, typecount);
1026
1027             for (tpi = 0; tpi < ntp; tpi++)
1028             {
1029                 for (tpj = tpi; tpj < ntp; tpj++)
1030                 {
1031                     tmpi = typecount[tpi];
1032                     tmpj = typecount[tpj];
1033                     if (tpi != tpj)
1034                     {
1035                         npair_ij = tmpi*tmpj;
1036                     }
1037                     else
1038                     {
1039                         npair_ij = tmpi*(tmpi - 1)/2;
1040                     }
1041                     if (bBHAM)
1042                     {
1043                         /* nbfp now includes the 6.0 derivative prefactor */
1044                         csix    += npair_ij*BHAMC(nbfp, ntp, tpi, tpj)/6.0;
1045                     }
1046                     else
1047                     {
1048                         /* nbfp now includes the 6.0/12.0 derivative prefactors */
1049                         csix    += npair_ij*   C6(nbfp, ntp, tpi, tpj)/6.0;
1050                         ctwelve += npair_ij*  C12(nbfp, ntp, tpi, tpj)/12.0;
1051                     }
1052                     npair += npair_ij;
1053                 }
1054             }
1055             sfree(typecount);
1056             /* Subtract the excluded pairs.
1057              * The main reason for substracting exclusions is that in some cases
1058              * some combinations might never occur and the parameters could have
1059              * any value. These unused values should not influence the dispersion
1060              * correction.
1061              */
1062             for (mb = 0; mb < mtop->nmolblock; mb++)
1063             {
1064                 nmol  = mtop->molblock[mb].nmol;
1065                 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
1066                 excl  = &mtop->moltype[mtop->molblock[mb].type].excls;
1067                 for (i = 0; (i < atoms->nr); i++)
1068                 {
1069                     if (q == 0)
1070                     {
1071                         tpi = atoms->atom[i].type;
1072                     }
1073                     else
1074                     {
1075                         tpi = atoms->atom[i].typeB;
1076                     }
1077                     j1  = excl->index[i];
1078                     j2  = excl->index[i+1];
1079                     for (j = j1; j < j2; j++)
1080                     {
1081                         k = excl->a[j];
1082                         if (k > i)
1083                         {
1084                             if (q == 0)
1085                             {
1086                                 tpj = atoms->atom[k].type;
1087                             }
1088                             else
1089                             {
1090                                 tpj = atoms->atom[k].typeB;
1091                             }
1092                             if (bBHAM)
1093                             {
1094                                 /* nbfp now includes the 6.0 derivative prefactor */
1095                                 csix -= nmol*BHAMC(nbfp, ntp, tpi, tpj)/6.0;
1096                             }
1097                             else
1098                             {
1099                                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1100                                 csix    -= nmol*C6 (nbfp, ntp, tpi, tpj)/6.0;
1101                                 ctwelve -= nmol*C12(nbfp, ntp, tpi, tpj)/12.0;
1102                             }
1103                             nexcl += nmol;
1104                         }
1105                     }
1106                 }
1107             }
1108         }
1109         else
1110         {
1111             /* Only correct for the interaction of the test particle
1112              * with the rest of the system.
1113              */
1114             atoms_tpi =
1115                 &mtop->moltype[mtop->molblock[mtop->nmolblock-1].type].atoms;
1116
1117             npair = 0;
1118             for (mb = 0; mb < mtop->nmolblock; mb++)
1119             {
1120                 nmol  = mtop->molblock[mb].nmol;
1121                 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
1122                 for (j = 0; j < atoms->nr; j++)
1123                 {
1124                     nmolc = nmol;
1125                     /* Remove the interaction of the test charge group
1126                      * with itself.
1127                      */
1128                     if (mb == mtop->nmolblock-1)
1129                     {
1130                         nmolc--;
1131
1132                         if (mb == 0 && nmol == 1)
1133                         {
1134                             gmx_fatal(FARGS, "Old format tpr with TPI, please generate a new tpr file");
1135                         }
1136                     }
1137                     if (q == 0)
1138                     {
1139                         tpj = atoms->atom[j].type;
1140                     }
1141                     else
1142                     {
1143                         tpj = atoms->atom[j].typeB;
1144                     }
1145                     for (i = 0; i < fr->n_tpi; i++)
1146                     {
1147                         if (q == 0)
1148                         {
1149                             tpi = atoms_tpi->atom[i].type;
1150                         }
1151                         else
1152                         {
1153                             tpi = atoms_tpi->atom[i].typeB;
1154                         }
1155                         if (bBHAM)
1156                         {
1157                             /* nbfp now includes the 6.0 derivative prefactor */
1158                             csix    += nmolc*BHAMC(nbfp, ntp, tpi, tpj)/6.0;
1159                         }
1160                         else
1161                         {
1162                             /* nbfp now includes the 6.0/12.0 derivative prefactors */
1163                             csix    += nmolc*C6 (nbfp, ntp, tpi, tpj)/6.0;
1164                             ctwelve += nmolc*C12(nbfp, ntp, tpi, tpj)/12.0;
1165                         }
1166                         npair += nmolc;
1167                     }
1168                 }
1169             }
1170         }
1171         if (npair - nexcl <= 0 && fplog)
1172         {
1173             fprintf(fplog, "\nWARNING: There are no atom pairs for dispersion correction\n\n");
1174             csix     = 0;
1175             ctwelve  = 0;
1176         }
1177         else
1178         {
1179             csix    /= npair - nexcl;
1180             ctwelve /= npair - nexcl;
1181         }
1182         if (debug)
1183         {
1184             fprintf(debug, "Counted %d exclusions\n", nexcl);
1185             fprintf(debug, "Average C6 parameter is: %10g\n", (double)csix);
1186             fprintf(debug, "Average C12 parameter is: %10g\n", (double)ctwelve);
1187         }
1188         fr->avcsix[q]    = csix;
1189         fr->avctwelve[q] = ctwelve;
1190     }
1191
1192     if (EVDW_PME(fr->vdwtype))
1193     {
1194         sfree(nbfp_comb);
1195     }
1196
1197     if (fplog != NULL)
1198     {
1199         if (fr->eDispCorr == edispcAllEner ||
1200             fr->eDispCorr == edispcAllEnerPres)
1201         {
1202             fprintf(fplog, "Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
1203                     fr->avcsix[0], fr->avctwelve[0]);
1204         }
1205         else
1206         {
1207             fprintf(fplog, "Long Range LJ corr.: <C6> %10.4e\n", fr->avcsix[0]);
1208         }
1209     }
1210 }
1211
1212
1213 static void set_bham_b_max(FILE *fplog, t_forcerec *fr,
1214                            const gmx_mtop_t *mtop)
1215 {
1216     const t_atoms *at1, *at2;
1217     int            mt1, mt2, i, j, tpi, tpj, ntypes;
1218     real           b, bmin;
1219     real          *nbfp;
1220
1221     if (fplog)
1222     {
1223         fprintf(fplog, "Determining largest Buckingham b parameter for table\n");
1224     }
1225     nbfp   = fr->nbfp;
1226     ntypes = fr->ntype;
1227
1228     bmin           = -1;
1229     fr->bham_b_max = 0;
1230     for (mt1 = 0; mt1 < mtop->nmoltype; mt1++)
1231     {
1232         at1 = &mtop->moltype[mt1].atoms;
1233         for (i = 0; (i < at1->nr); i++)
1234         {
1235             tpi = at1->atom[i].type;
1236             if (tpi >= ntypes)
1237             {
1238                 gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes);
1239             }
1240
1241             for (mt2 = mt1; mt2 < mtop->nmoltype; mt2++)
1242             {
1243                 at2 = &mtop->moltype[mt2].atoms;
1244                 for (j = 0; (j < at2->nr); j++)
1245                 {
1246                     tpj = at2->atom[j].type;
1247                     if (tpj >= ntypes)
1248                     {
1249                         gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes);
1250                     }
1251                     b = BHAMB(nbfp, ntypes, tpi, tpj);
1252                     if (b > fr->bham_b_max)
1253                     {
1254                         fr->bham_b_max = b;
1255                     }
1256                     if ((b < bmin) || (bmin == -1))
1257                     {
1258                         bmin = b;
1259                     }
1260                 }
1261             }
1262         }
1263     }
1264     if (fplog)
1265     {
1266         fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n",
1267                 bmin, fr->bham_b_max);
1268     }
1269 }
1270
1271 static void make_nbf_tables(FILE *fp, const output_env_t oenv,
1272                             t_forcerec *fr, real rtab,
1273                             const t_commrec *cr,
1274                             const char *tabfn, char *eg1, char *eg2,
1275                             t_nblists *nbl)
1276 {
1277     char buf[STRLEN];
1278     int  i, j;
1279
1280     if (tabfn == NULL)
1281     {
1282         if (debug)
1283         {
1284             fprintf(debug, "No table file name passed, can not read table, can not do non-bonded interactions\n");
1285         }
1286         return;
1287     }
1288
1289     sprintf(buf, "%s", tabfn);
1290     if (eg1 && eg2)
1291     {
1292         /* Append the two energy group names */
1293         sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "_%s_%s.%s",
1294                 eg1, eg2, ftp2ext(efXVG));
1295     }
1296     nbl->table_elec_vdw = make_tables(fp, oenv, fr, MASTER(cr), buf, rtab, 0);
1297     /* Copy the contents of the table to separate coulomb and LJ tables too,
1298      * to improve cache performance.
1299      */
1300     /* For performance reasons we want
1301      * the table data to be aligned to 16-byte. The pointers could be freed
1302      * but currently aren't.
1303      */
1304     nbl->table_elec.interaction   = GMX_TABLE_INTERACTION_ELEC;
1305     nbl->table_elec.format        = nbl->table_elec_vdw.format;
1306     nbl->table_elec.r             = nbl->table_elec_vdw.r;
1307     nbl->table_elec.n             = nbl->table_elec_vdw.n;
1308     nbl->table_elec.scale         = nbl->table_elec_vdw.scale;
1309     nbl->table_elec.scale_exp     = nbl->table_elec_vdw.scale_exp;
1310     nbl->table_elec.formatsize    = nbl->table_elec_vdw.formatsize;
1311     nbl->table_elec.ninteractions = 1;
1312     nbl->table_elec.stride        = nbl->table_elec.formatsize * nbl->table_elec.ninteractions;
1313     snew_aligned(nbl->table_elec.data, nbl->table_elec.stride*(nbl->table_elec.n+1), 32);
1314
1315     nbl->table_vdw.interaction   = GMX_TABLE_INTERACTION_VDWREP_VDWDISP;
1316     nbl->table_vdw.format        = nbl->table_elec_vdw.format;
1317     nbl->table_vdw.r             = nbl->table_elec_vdw.r;
1318     nbl->table_vdw.n             = nbl->table_elec_vdw.n;
1319     nbl->table_vdw.scale         = nbl->table_elec_vdw.scale;
1320     nbl->table_vdw.scale_exp     = nbl->table_elec_vdw.scale_exp;
1321     nbl->table_vdw.formatsize    = nbl->table_elec_vdw.formatsize;
1322     nbl->table_vdw.ninteractions = 2;
1323     nbl->table_vdw.stride        = nbl->table_vdw.formatsize * nbl->table_vdw.ninteractions;
1324     snew_aligned(nbl->table_vdw.data, nbl->table_vdw.stride*(nbl->table_vdw.n+1), 32);
1325
1326     for (i = 0; i <= nbl->table_elec_vdw.n; i++)
1327     {
1328         for (j = 0; j < 4; j++)
1329         {
1330             nbl->table_elec.data[4*i+j] = nbl->table_elec_vdw.data[12*i+j];
1331         }
1332         for (j = 0; j < 8; j++)
1333         {
1334             nbl->table_vdw.data[8*i+j] = nbl->table_elec_vdw.data[12*i+4+j];
1335         }
1336     }
1337 }
1338
1339 static void count_tables(int ftype1, int ftype2, const gmx_mtop_t *mtop,
1340                          int *ncount, int **count)
1341 {
1342     const gmx_moltype_t *molt;
1343     const t_ilist       *il;
1344     int                  mt, ftype, stride, i, j, tabnr;
1345
1346     for (mt = 0; mt < mtop->nmoltype; mt++)
1347     {
1348         molt = &mtop->moltype[mt];
1349         for (ftype = 0; ftype < F_NRE; ftype++)
1350         {
1351             if (ftype == ftype1 || ftype == ftype2)
1352             {
1353                 il     = &molt->ilist[ftype];
1354                 stride = 1 + NRAL(ftype);
1355                 for (i = 0; i < il->nr; i += stride)
1356                 {
1357                     tabnr = mtop->ffparams.iparams[il->iatoms[i]].tab.table;
1358                     if (tabnr < 0)
1359                     {
1360                         gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr);
1361                     }
1362                     if (tabnr >= *ncount)
1363                     {
1364                         srenew(*count, tabnr+1);
1365                         for (j = *ncount; j < tabnr+1; j++)
1366                         {
1367                             (*count)[j] = 0;
1368                         }
1369                         *ncount = tabnr+1;
1370                     }
1371                     (*count)[tabnr]++;
1372                 }
1373             }
1374         }
1375     }
1376 }
1377
1378 static bondedtable_t *make_bonded_tables(FILE *fplog,
1379                                          int ftype1, int ftype2,
1380                                          const gmx_mtop_t *mtop,
1381                                          const char *basefn, const char *tabext)
1382 {
1383     int            i, ncount, *count;
1384     char           tabfn[STRLEN];
1385     bondedtable_t *tab;
1386
1387     tab = NULL;
1388
1389     ncount = 0;
1390     count  = NULL;
1391     count_tables(ftype1, ftype2, mtop, &ncount, &count);
1392
1393     if (ncount > 0)
1394     {
1395         snew(tab, ncount);
1396         for (i = 0; i < ncount; i++)
1397         {
1398             if (count[i] > 0)
1399             {
1400                 sprintf(tabfn, "%s", basefn);
1401                 sprintf(tabfn + strlen(basefn) - strlen(ftp2ext(efXVG)) - 1, "_%s%d.%s",
1402                         tabext, i, ftp2ext(efXVG));
1403                 tab[i] = make_bonded_table(fplog, tabfn, NRAL(ftype1)-2);
1404             }
1405         }
1406         sfree(count);
1407     }
1408
1409     return tab;
1410 }
1411
1412 void forcerec_set_ranges(t_forcerec *fr,
1413                          int ncg_home, int ncg_force,
1414                          int natoms_force,
1415                          int natoms_force_constr, int natoms_f_novirsum)
1416 {
1417     fr->cg0 = 0;
1418     fr->hcg = ncg_home;
1419
1420     /* fr->ncg_force is unused in the standard code,
1421      * but it can be useful for modified code dealing with charge groups.
1422      */
1423     fr->ncg_force           = ncg_force;
1424     fr->natoms_force        = natoms_force;
1425     fr->natoms_force_constr = natoms_force_constr;
1426
1427     if (fr->natoms_force_constr > fr->nalloc_force)
1428     {
1429         fr->nalloc_force = over_alloc_dd(fr->natoms_force_constr);
1430
1431         if (fr->bTwinRange)
1432         {
1433             srenew(fr->f_twin, fr->nalloc_force);
1434         }
1435     }
1436
1437     if (fr->bF_NoVirSum)
1438     {
1439         fr->f_novirsum_n = natoms_f_novirsum;
1440         if (fr->f_novirsum_n > fr->f_novirsum_nalloc)
1441         {
1442             fr->f_novirsum_nalloc = over_alloc_dd(fr->f_novirsum_n);
1443             srenew(fr->f_novirsum_alloc, fr->f_novirsum_nalloc);
1444         }
1445     }
1446     else
1447     {
1448         fr->f_novirsum_n = 0;
1449     }
1450 }
1451
1452 static real cutoff_inf(real cutoff)
1453 {
1454     if (cutoff == 0)
1455     {
1456         cutoff = GMX_CUTOFF_INF;
1457     }
1458
1459     return cutoff;
1460 }
1461
1462 static void make_adress_tf_tables(FILE *fp, const output_env_t oenv,
1463                                   t_forcerec *fr, const t_inputrec *ir,
1464                                   const char *tabfn, const gmx_mtop_t *mtop,
1465                                   matrix     box)
1466 {
1467     char buf[STRLEN];
1468     int  i, j;
1469
1470     if (tabfn == NULL)
1471     {
1472         gmx_fatal(FARGS, "No thermoforce table file given. Use -tabletf to specify a file\n");
1473         return;
1474     }
1475
1476     snew(fr->atf_tabs, ir->adress->n_tf_grps);
1477
1478     sprintf(buf, "%s", tabfn);
1479     for (i = 0; i < ir->adress->n_tf_grps; i++)
1480     {
1481         j = ir->adress->tf_table_index[i]; /* get energy group index */
1482         sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "tf_%s.%s",
1483                 *(mtop->groups.grpname[mtop->groups.grps[egcENER].nm_ind[j]]), ftp2ext(efXVG));
1484         if (fp)
1485         {
1486             fprintf(fp, "loading tf table for energygrp index %d from %s\n", ir->adress->tf_table_index[i], buf);
1487         }
1488         fr->atf_tabs[i] = make_atf_table(fp, oenv, fr, buf, box);
1489     }
1490
1491 }
1492
1493 gmx_bool can_use_allvsall(const t_inputrec *ir, gmx_bool bPrintNote, t_commrec *cr, FILE *fp)
1494 {
1495     gmx_bool bAllvsAll;
1496
1497     bAllvsAll =
1498         (
1499             ir->rlist == 0            &&
1500             ir->rcoulomb == 0         &&
1501             ir->rvdw == 0             &&
1502             ir->ePBC == epbcNONE      &&
1503             ir->vdwtype == evdwCUT    &&
1504             ir->coulombtype == eelCUT &&
1505             ir->efep == efepNO        &&
1506             (ir->implicit_solvent == eisNO ||
1507              (ir->implicit_solvent == eisGBSA && (ir->gb_algorithm == egbSTILL ||
1508                                                   ir->gb_algorithm == egbHCT   ||
1509                                                   ir->gb_algorithm == egbOBC))) &&
1510             getenv("GMX_NO_ALLVSALL") == NULL
1511         );
1512
1513     if (bAllvsAll && ir->opts.ngener > 1)
1514     {
1515         const char *note = "NOTE: Can not use all-vs-all force loops, because there are multiple energy monitor groups; you might get significantly higher performance when using only a single energy monitor group.\n";
1516
1517         if (bPrintNote)
1518         {
1519             if (MASTER(cr))
1520             {
1521                 fprintf(stderr, "\n%s\n", note);
1522             }
1523             if (fp != NULL)
1524             {
1525                 fprintf(fp, "\n%s\n", note);
1526             }
1527         }
1528         bAllvsAll = FALSE;
1529     }
1530
1531     if (bAllvsAll && fp && MASTER(cr))
1532     {
1533         fprintf(fp, "\nUsing SIMD all-vs-all kernels.\n\n");
1534     }
1535
1536     return bAllvsAll;
1537 }
1538
1539
1540 static void init_forcerec_f_threads(t_forcerec *fr, int nenergrp)
1541 {
1542     int t, i;
1543
1544     /* These thread local data structures are used for bondeds only */
1545     fr->nthreads = gmx_omp_nthreads_get(emntBonded);
1546
1547     if (fr->nthreads > 1)
1548     {
1549         snew(fr->f_t, fr->nthreads);
1550         /* Thread 0 uses the global force and energy arrays */
1551         for (t = 1; t < fr->nthreads; t++)
1552         {
1553             fr->f_t[t].f        = NULL;
1554             fr->f_t[t].f_nalloc = 0;
1555             snew(fr->f_t[t].fshift, SHIFTS);
1556             fr->f_t[t].grpp.nener = nenergrp*nenergrp;
1557             for (i = 0; i < egNR; i++)
1558             {
1559                 snew(fr->f_t[t].grpp.ener[i], fr->f_t[t].grpp.nener);
1560             }
1561         }
1562     }
1563 }
1564
1565
1566 gmx_bool nbnxn_acceleration_supported(FILE             *fplog,
1567                                       const t_commrec  *cr,
1568                                       const t_inputrec *ir,
1569                                       gmx_bool          bGPU)
1570 {
1571     if (!bGPU && (ir->vdwtype == evdwPME && ir->ljpme_combination_rule == eljpmeLB))
1572     {
1573         md_print_warn(cr, fplog, "LJ-PME with Lorentz-Berthelot is not supported with %s, falling back to %s\n",
1574                       bGPU ? "GPUs" : "SIMD kernels",
1575                       bGPU ? "CPU only" : "plain-C kernels");
1576         return FALSE;
1577     }
1578
1579     return TRUE;
1580 }
1581
1582
1583 static void pick_nbnxn_kernel_cpu(const t_inputrec gmx_unused *ir,
1584                                   int                         *kernel_type,
1585                                   int                         *ewald_excl)
1586 {
1587     *kernel_type = nbnxnk4x4_PlainC;
1588     *ewald_excl  = ewaldexclTable;
1589
1590 #ifdef GMX_NBNXN_SIMD
1591     {
1592 #ifdef GMX_NBNXN_SIMD_4XN
1593         *kernel_type = nbnxnk4xN_SIMD_4xN;
1594 #endif
1595 #ifdef GMX_NBNXN_SIMD_2XNN
1596         *kernel_type = nbnxnk4xN_SIMD_2xNN;
1597 #endif
1598
1599 #if defined GMX_NBNXN_SIMD_2XNN && defined GMX_NBNXN_SIMD_4XN
1600         /* We need to choose if we want 2x(N+N) or 4xN kernels.
1601          * Currently this is based on the SIMD acceleration choice,
1602          * but it might be better to decide this at runtime based on CPU.
1603          *
1604          * 4xN calculates more (zero) interactions, but has less pair-search
1605          * work and much better kernel instruction scheduling.
1606          *
1607          * Up till now we have only seen that on Intel Sandy/Ivy Bridge,
1608          * which doesn't have FMA, both the analytical and tabulated Ewald
1609          * kernels have similar pair rates for 4x8 and 2x(4+4), so we choose
1610          * 2x(4+4) because it results in significantly fewer pairs.
1611          * For RF, the raw pair rate of the 4x8 kernel is higher than 2x(4+4),
1612          * 10% with HT, 50% without HT. As we currently don't detect the actual
1613          * use of HT, use 4x8 to avoid a potential performance hit.
1614          * On Intel Haswell 4x8 is always faster.
1615          */
1616         *kernel_type = nbnxnk4xN_SIMD_4xN;
1617
1618 #ifndef GMX_SIMD_HAVE_FMA
1619         if (EEL_PME_EWALD(ir->coulombtype) ||
1620             EVDW_PME(ir->vdwtype))
1621         {
1622             /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
1623              * There are enough instructions to make 2x(4+4) efficient.
1624              */
1625             *kernel_type = nbnxnk4xN_SIMD_2xNN;
1626         }
1627 #endif
1628 #endif  /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
1629
1630
1631         if (getenv("GMX_NBNXN_SIMD_4XN") != NULL)
1632         {
1633 #ifdef GMX_NBNXN_SIMD_4XN
1634             *kernel_type = nbnxnk4xN_SIMD_4xN;
1635 #else
1636             gmx_fatal(FARGS, "SIMD 4xN kernels requested, but Gromacs has been compiled without support for these kernels");
1637 #endif
1638         }
1639         if (getenv("GMX_NBNXN_SIMD_2XNN") != NULL)
1640         {
1641 #ifdef GMX_NBNXN_SIMD_2XNN
1642             *kernel_type = nbnxnk4xN_SIMD_2xNN;
1643 #else
1644             gmx_fatal(FARGS, "SIMD 2x(N+N) kernels requested, but Gromacs has been compiled without support for these kernels");
1645 #endif
1646         }
1647
1648         /* Analytical Ewald exclusion correction is only an option in
1649          * the SIMD kernel.
1650          * Since table lookup's don't parallelize with SIMD, analytical
1651          * will probably always be faster for a SIMD width of 8 or more.
1652          * With FMA analytical is sometimes faster for a width if 4 as well.
1653          * On BlueGene/Q, this is faster regardless of precision.
1654          * In single precision, this is faster on Bulldozer.
1655          */
1656 #if GMX_SIMD_REAL_WIDTH >= 8 || \
1657         (GMX_SIMD_REAL_WIDTH >= 4 && defined GMX_SIMD_HAVE_FMA && !defined GMX_DOUBLE) || \
1658         defined GMX_SIMD_IBM_QPX
1659         *ewald_excl = ewaldexclAnalytical;
1660 #endif
1661         if (getenv("GMX_NBNXN_EWALD_TABLE") != NULL)
1662         {
1663             *ewald_excl = ewaldexclTable;
1664         }
1665         if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != NULL)
1666         {
1667             *ewald_excl = ewaldexclAnalytical;
1668         }
1669
1670     }
1671 #endif /* GMX_NBNXN_SIMD */
1672 }
1673
1674
1675 const char *lookup_nbnxn_kernel_name(int kernel_type)
1676 {
1677     const char *returnvalue = NULL;
1678     switch (kernel_type)
1679     {
1680         case nbnxnkNotSet:
1681             returnvalue = "not set";
1682             break;
1683         case nbnxnk4x4_PlainC:
1684             returnvalue = "plain C";
1685             break;
1686         case nbnxnk4xN_SIMD_4xN:
1687         case nbnxnk4xN_SIMD_2xNN:
1688 #ifdef GMX_NBNXN_SIMD
1689 #if defined GMX_SIMD_X86_SSE2
1690             returnvalue = "SSE2";
1691 #elif defined GMX_SIMD_X86_SSE4_1
1692             returnvalue = "SSE4.1";
1693 #elif defined GMX_SIMD_X86_AVX_128_FMA
1694             returnvalue = "AVX_128_FMA";
1695 #elif defined GMX_SIMD_X86_AVX_256
1696             returnvalue = "AVX_256";
1697 #elif defined GMX_SIMD_X86_AVX2_256
1698             returnvalue = "AVX2_256";
1699 #else
1700             returnvalue = "SIMD";
1701 #endif
1702 #else  /* GMX_NBNXN_SIMD */
1703             returnvalue = "not available";
1704 #endif /* GMX_NBNXN_SIMD */
1705             break;
1706         case nbnxnk8x8x8_CUDA: returnvalue   = "CUDA"; break;
1707         case nbnxnk8x8x8_PlainC: returnvalue = "plain C"; break;
1708
1709         case nbnxnkNR:
1710         default:
1711             gmx_fatal(FARGS, "Illegal kernel type selected");
1712             returnvalue = NULL;
1713             break;
1714     }
1715     return returnvalue;
1716 };
1717
1718 static void pick_nbnxn_kernel(FILE                *fp,
1719                               const t_commrec     *cr,
1720                               gmx_bool             use_simd_kernels,
1721                               gmx_bool             bUseGPU,
1722                               gmx_bool             bEmulateGPU,
1723                               const t_inputrec    *ir,
1724                               int                 *kernel_type,
1725                               int                 *ewald_excl,
1726                               gmx_bool             bDoNonbonded)
1727 {
1728     assert(kernel_type);
1729
1730     *kernel_type = nbnxnkNotSet;
1731     *ewald_excl  = ewaldexclTable;
1732
1733     if (bEmulateGPU)
1734     {
1735         *kernel_type = nbnxnk8x8x8_PlainC;
1736
1737         if (bDoNonbonded)
1738         {
1739             md_print_warn(cr, fp, "Emulating a GPU run on the CPU (slow)");
1740         }
1741     }
1742     else if (bUseGPU)
1743     {
1744         *kernel_type = nbnxnk8x8x8_CUDA;
1745     }
1746
1747     if (*kernel_type == nbnxnkNotSet)
1748     {
1749         /* LJ PME with LB combination rule does 7 mesh operations.
1750          * This so slow that we don't compile SIMD non-bonded kernels for that.
1751          */
1752         if (use_simd_kernels &&
1753             nbnxn_acceleration_supported(fp, cr, ir, FALSE))
1754         {
1755             pick_nbnxn_kernel_cpu(ir, kernel_type, ewald_excl);
1756         }
1757         else
1758         {
1759             *kernel_type = nbnxnk4x4_PlainC;
1760         }
1761     }
1762
1763     if (bDoNonbonded && fp != NULL)
1764     {
1765         fprintf(fp, "\nUsing %s %dx%d non-bonded kernels\n\n",
1766                 lookup_nbnxn_kernel_name(*kernel_type),
1767                 nbnxn_kernel_pairlist_simple(*kernel_type) ? NBNXN_CPU_CLUSTER_I_SIZE : NBNXN_GPU_CLUSTER_SIZE,
1768                 nbnxn_kernel_to_cj_size(*kernel_type));
1769
1770         if (nbnxnk4x4_PlainC == *kernel_type ||
1771             nbnxnk8x8x8_PlainC == *kernel_type)
1772         {
1773             md_print_warn(cr, fp,
1774                           "WARNING: Using the slow %s kernels. This should\n"
1775                           "not happen during routine usage on supported platforms.\n\n",
1776                           lookup_nbnxn_kernel_name(*kernel_type));
1777         }
1778     }
1779 }
1780
1781 static void pick_nbnxn_resources(FILE                *fp,
1782                                  const t_commrec     *cr,
1783                                  const gmx_hw_info_t *hwinfo,
1784                                  gmx_bool             bDoNonbonded,
1785                                  gmx_bool            *bUseGPU,
1786                                  gmx_bool            *bEmulateGPU,
1787                                  const gmx_gpu_opt_t *gpu_opt)
1788 {
1789     gmx_bool bEmulateGPUEnvVarSet;
1790     char     gpu_err_str[STRLEN];
1791
1792     *bUseGPU = FALSE;
1793
1794     bEmulateGPUEnvVarSet = (getenv("GMX_EMULATE_GPU") != NULL);
1795
1796     /* Run GPU emulation mode if GMX_EMULATE_GPU is defined. Because
1797      * GPUs (currently) only handle non-bonded calculations, we will
1798      * automatically switch to emulation if non-bonded calculations are
1799      * turned off via GMX_NO_NONBONDED - this is the simple and elegant
1800      * way to turn off GPU initialization, data movement, and cleanup.
1801      *
1802      * GPU emulation can be useful to assess the performance one can expect by
1803      * adding GPU(s) to the machine. The conditional below allows this even
1804      * if mdrun is compiled without GPU acceleration support.
1805      * Note that you should freezing the system as otherwise it will explode.
1806      */
1807     *bEmulateGPU = (bEmulateGPUEnvVarSet ||
1808                     (!bDoNonbonded &&
1809                      gpu_opt->ncuda_dev_use > 0));
1810
1811     /* Enable GPU mode when GPUs are available or no GPU emulation is requested.
1812      */
1813     if (gpu_opt->ncuda_dev_use > 0 && !(*bEmulateGPU))
1814     {
1815         /* Each PP node will use the intra-node id-th device from the
1816          * list of detected/selected GPUs. */
1817         if (!init_gpu(fp, cr->rank_pp_intranode, gpu_err_str,
1818                       &hwinfo->gpu_info, gpu_opt))
1819         {
1820             /* At this point the init should never fail as we made sure that
1821              * we have all the GPUs we need. If it still does, we'll bail. */
1822             gmx_fatal(FARGS, "On rank %d failed to initialize GPU #%d: %s",
1823                       cr->nodeid,
1824                       get_gpu_device_id(&hwinfo->gpu_info, gpu_opt,
1825                                         cr->rank_pp_intranode),
1826                       gpu_err_str);
1827         }
1828
1829         /* Here we actually turn on hardware GPU acceleration */
1830         *bUseGPU = TRUE;
1831     }
1832 }
1833
1834 gmx_bool uses_simple_tables(int                 cutoff_scheme,
1835                             nonbonded_verlet_t *nbv,
1836                             int                 group)
1837 {
1838     gmx_bool bUsesSimpleTables = TRUE;
1839     int      grp_index;
1840
1841     switch (cutoff_scheme)
1842     {
1843         case ecutsGROUP:
1844             bUsesSimpleTables = TRUE;
1845             break;
1846         case ecutsVERLET:
1847             assert(NULL != nbv && NULL != nbv->grp);
1848             grp_index         = (group < 0) ? 0 : (nbv->ngrp - 1);
1849             bUsesSimpleTables = nbnxn_kernel_pairlist_simple(nbv->grp[grp_index].kernel_type);
1850             break;
1851         default:
1852             gmx_incons("unimplemented");
1853     }
1854     return bUsesSimpleTables;
1855 }
1856
1857 static void init_ewald_f_table(interaction_const_t *ic,
1858                                gmx_bool             bUsesSimpleTables,
1859                                real                 rtab)
1860 {
1861     real maxr;
1862
1863     if (bUsesSimpleTables)
1864     {
1865         /* Get the Ewald table spacing based on Coulomb and/or LJ
1866          * Ewald coefficients and rtol.
1867          */
1868         ic->tabq_scale = ewald_spline3_table_scale(ic);
1869
1870         maxr           = (rtab > ic->rcoulomb) ? rtab : ic->rcoulomb;
1871         ic->tabq_size  = (int)(maxr*ic->tabq_scale) + 2;
1872     }
1873     else
1874     {
1875         ic->tabq_size = GPU_EWALD_COULOMB_FORCE_TABLE_SIZE;
1876         /* Subtract 2 iso 1 to avoid access out of range due to rounding */
1877         ic->tabq_scale = (ic->tabq_size - 2)/ic->rcoulomb;
1878     }
1879
1880     sfree_aligned(ic->tabq_coul_FDV0);
1881     sfree_aligned(ic->tabq_coul_F);
1882     sfree_aligned(ic->tabq_coul_V);
1883
1884     sfree_aligned(ic->tabq_vdw_FDV0);
1885     sfree_aligned(ic->tabq_vdw_F);
1886     sfree_aligned(ic->tabq_vdw_V);
1887
1888     if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype))
1889     {
1890         /* Create the original table data in FDV0 */
1891         snew_aligned(ic->tabq_coul_FDV0, ic->tabq_size*4, 32);
1892         snew_aligned(ic->tabq_coul_F, ic->tabq_size, 32);
1893         snew_aligned(ic->tabq_coul_V, ic->tabq_size, 32);
1894         table_spline3_fill_ewald_lr(ic->tabq_coul_F, ic->tabq_coul_V, ic->tabq_coul_FDV0,
1895                                     ic->tabq_size, 1/ic->tabq_scale, ic->ewaldcoeff_q, v_q_ewald_lr);
1896     }
1897
1898     if (EVDW_PME(ic->vdwtype))
1899     {
1900         snew_aligned(ic->tabq_vdw_FDV0, ic->tabq_size*4, 32);
1901         snew_aligned(ic->tabq_vdw_F, ic->tabq_size, 32);
1902         snew_aligned(ic->tabq_vdw_V, ic->tabq_size, 32);
1903         table_spline3_fill_ewald_lr(ic->tabq_vdw_F, ic->tabq_vdw_V, ic->tabq_vdw_FDV0,
1904                                     ic->tabq_size, 1/ic->tabq_scale, ic->ewaldcoeff_lj, v_lj_ewald_lr);
1905     }
1906 }
1907
1908 void init_interaction_const_tables(FILE                *fp,
1909                                    interaction_const_t *ic,
1910                                    gmx_bool             bUsesSimpleTables,
1911                                    real                 rtab)
1912 {
1913     if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype) || EVDW_PME(ic->vdwtype))
1914     {
1915         init_ewald_f_table(ic, bUsesSimpleTables, rtab);
1916
1917         if (fp != NULL)
1918         {
1919             fprintf(fp, "Initialized non-bonded Ewald correction tables, spacing: %.2e size: %d\n\n",
1920                     1/ic->tabq_scale, ic->tabq_size);
1921         }
1922     }
1923 }
1924
1925 static void clear_force_switch_constants(shift_consts_t *sc)
1926 {
1927     sc->c2   = 0;
1928     sc->c3   = 0;
1929     sc->cpot = 0;
1930 }
1931
1932 static void force_switch_constants(real p,
1933                                    real rsw, real rc,
1934                                    shift_consts_t *sc)
1935 {
1936     /* Here we determine the coefficient for shifting the force to zero
1937      * between distance rsw and the cut-off rc.
1938      * For a potential of r^-p, we have force p*r^-(p+1).
1939      * But to save flops we absorb p in the coefficient.
1940      * Thus we get:
1941      * force/p   = r^-(p+1) + c2*r^2 + c3*r^3
1942      * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot
1943      */
1944     sc->c2   =  ((p + 1)*rsw - (p + 4)*rc)/(pow(rc, p + 2)*pow(rc - rsw, 2));
1945     sc->c3   = -((p + 1)*rsw - (p + 3)*rc)/(pow(rc, p + 2)*pow(rc - rsw, 3));
1946     sc->cpot = -pow(rc, -p) + p*sc->c2/3*pow(rc - rsw, 3) + p*sc->c3/4*pow(rc - rsw, 4);
1947 }
1948
1949 static void potential_switch_constants(real rsw, real rc,
1950                                        switch_consts_t *sc)
1951 {
1952     /* The switch function is 1 at rsw and 0 at rc.
1953      * The derivative and second derivate are zero at both ends.
1954      * rsw        = max(r - r_switch, 0)
1955      * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
1956      * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
1957      * force      = force*dsw - potential*sw
1958      * potential *= sw
1959      */
1960     sc->c3 = -10*pow(rc - rsw, -3);
1961     sc->c4 =  15*pow(rc - rsw, -4);
1962     sc->c5 =  -6*pow(rc - rsw, -5);
1963 }
1964
1965 static void
1966 init_interaction_const(FILE                       *fp,
1967                        const t_commrec gmx_unused *cr,
1968                        interaction_const_t       **interaction_const,
1969                        const t_forcerec           *fr,
1970                        real                        rtab)
1971 {
1972     interaction_const_t *ic;
1973     gmx_bool             bUsesSimpleTables = TRUE;
1974     const real           minusSix          = -6.0;
1975     const real           minusTwelve       = -12.0;
1976
1977     snew(ic, 1);
1978
1979     /* Just allocate something so we can free it */
1980     snew_aligned(ic->tabq_coul_FDV0, 16, 32);
1981     snew_aligned(ic->tabq_coul_F, 16, 32);
1982     snew_aligned(ic->tabq_coul_V, 16, 32);
1983
1984     ic->rlist           = fr->rlist;
1985     ic->rlistlong       = fr->rlistlong;
1986
1987     /* Lennard-Jones */
1988     ic->vdwtype         = fr->vdwtype;
1989     ic->vdw_modifier    = fr->vdw_modifier;
1990     ic->rvdw            = fr->rvdw;
1991     ic->rvdw_switch     = fr->rvdw_switch;
1992     ic->ewaldcoeff_lj   = fr->ewaldcoeff_lj;
1993     ic->ljpme_comb_rule = fr->ljpme_combination_rule;
1994     ic->sh_lj_ewald     = 0;
1995     clear_force_switch_constants(&ic->dispersion_shift);
1996     clear_force_switch_constants(&ic->repulsion_shift);
1997
1998     switch (ic->vdw_modifier)
1999     {
2000         case eintmodPOTSHIFT:
2001             /* Only shift the potential, don't touch the force */
2002             ic->dispersion_shift.cpot = -pow(ic->rvdw, minusSix);
2003             ic->repulsion_shift.cpot  = -pow(ic->rvdw, minusTwelve);
2004             if (EVDW_PME(ic->vdwtype))
2005             {
2006                 real crc2;
2007
2008                 crc2            = sqr(ic->ewaldcoeff_lj*ic->rvdw);
2009                 ic->sh_lj_ewald = (exp(-crc2)*(1 + crc2 + 0.5*crc2*crc2) - 1)*pow(ic->rvdw, minusSix);
2010             }
2011             break;
2012         case eintmodFORCESWITCH:
2013             /* Switch the force, switch and shift the potential */
2014             force_switch_constants(6.0, ic->rvdw_switch, ic->rvdw,
2015                                    &ic->dispersion_shift);
2016             force_switch_constants(12.0, ic->rvdw_switch, ic->rvdw,
2017                                    &ic->repulsion_shift);
2018             break;
2019         case eintmodPOTSWITCH:
2020             /* Switch the potential and force */
2021             potential_switch_constants(ic->rvdw_switch, ic->rvdw,
2022                                        &ic->vdw_switch);
2023             break;
2024         case eintmodNONE:
2025         case eintmodEXACTCUTOFF:
2026             /* Nothing to do here */
2027             break;
2028         default:
2029             gmx_incons("unimplemented potential modifier");
2030     }
2031
2032     ic->sh_invrc6 = -ic->dispersion_shift.cpot;
2033
2034     /* Electrostatics */
2035     ic->eeltype          = fr->eeltype;
2036     ic->coulomb_modifier = fr->coulomb_modifier;
2037     ic->rcoulomb         = fr->rcoulomb;
2038     ic->epsilon_r        = fr->epsilon_r;
2039     ic->epsfac           = fr->epsfac;
2040     ic->ewaldcoeff_q     = fr->ewaldcoeff_q;
2041
2042     if (fr->coulomb_modifier == eintmodPOTSHIFT)
2043     {
2044         ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
2045     }
2046     else
2047     {
2048         ic->sh_ewald = 0;
2049     }
2050
2051     /* Reaction-field */
2052     if (EEL_RF(ic->eeltype))
2053     {
2054         ic->epsilon_rf = fr->epsilon_rf;
2055         ic->k_rf       = fr->k_rf;
2056         ic->c_rf       = fr->c_rf;
2057     }
2058     else
2059     {
2060         /* For plain cut-off we might use the reaction-field kernels */
2061         ic->epsilon_rf = ic->epsilon_r;
2062         ic->k_rf       = 0;
2063         if (fr->coulomb_modifier == eintmodPOTSHIFT)
2064         {
2065             ic->c_rf   = 1/ic->rcoulomb;
2066         }
2067         else
2068         {
2069             ic->c_rf   = 0;
2070         }
2071     }
2072
2073     if (fp != NULL)
2074     {
2075         real dispersion_shift;
2076
2077         dispersion_shift = ic->dispersion_shift.cpot;
2078         if (EVDW_PME(ic->vdwtype))
2079         {
2080             dispersion_shift -= ic->sh_lj_ewald;
2081         }
2082         fprintf(fp, "Potential shift: LJ r^-12: %.3e r^-6: %.3e",
2083                 ic->repulsion_shift.cpot, dispersion_shift);
2084
2085         if (ic->eeltype == eelCUT)
2086         {
2087             fprintf(fp, ", Coulomb %.e", -ic->c_rf);
2088         }
2089         else if (EEL_PME(ic->eeltype))
2090         {
2091             fprintf(fp, ", Ewald %.3e", -ic->sh_ewald);
2092         }
2093         fprintf(fp, "\n");
2094     }
2095
2096     *interaction_const = ic;
2097
2098     if (fr->nbv != NULL && fr->nbv->bUseGPU)
2099     {
2100         nbnxn_cuda_init_const(fr->nbv->cu_nbv, ic, fr->nbv->grp);
2101
2102         /* With tMPI + GPUs some ranks may be sharing GPU(s) and therefore
2103          * also sharing texture references. To keep the code simple, we don't
2104          * treat texture references as shared resources, but this means that
2105          * the coulomb_tab and nbfp texture refs will get updated by multiple threads.
2106          * Hence, to ensure that the non-bonded kernels don't start before all
2107          * texture binding operations are finished, we need to wait for all ranks
2108          * to arrive here before continuing.
2109          *
2110          * Note that we could omit this barrier if GPUs are not shared (or
2111          * texture objects are used), but as this is initialization code, there
2112          * is not point in complicating things.
2113          */
2114 #ifdef GMX_THREAD_MPI
2115         if (PAR(cr))
2116         {
2117             gmx_barrier(cr);
2118         }
2119 #endif  /* GMX_THREAD_MPI */
2120     }
2121
2122     bUsesSimpleTables = uses_simple_tables(fr->cutoff_scheme, fr->nbv, -1);
2123     init_interaction_const_tables(fp, ic, bUsesSimpleTables, rtab);
2124 }
2125
2126 static void init_nb_verlet(FILE                *fp,
2127                            nonbonded_verlet_t **nb_verlet,
2128                            gmx_bool             bFEP_NonBonded,
2129                            const t_inputrec    *ir,
2130                            const t_forcerec    *fr,
2131                            const t_commrec     *cr,
2132                            const char          *nbpu_opt)
2133 {
2134     nonbonded_verlet_t *nbv;
2135     int                 i;
2136     char               *env;
2137     gmx_bool            bEmulateGPU, bHybridGPURun = FALSE;
2138
2139     nbnxn_alloc_t      *nb_alloc;
2140     nbnxn_free_t       *nb_free;
2141
2142     snew(nbv, 1);
2143
2144     pick_nbnxn_resources(fp, cr, fr->hwinfo,
2145                          fr->bNonbonded,
2146                          &nbv->bUseGPU,
2147                          &bEmulateGPU,
2148                          fr->gpu_opt);
2149
2150     nbv->nbs = NULL;
2151
2152     nbv->ngrp = (DOMAINDECOMP(cr) ? 2 : 1);
2153     for (i = 0; i < nbv->ngrp; i++)
2154     {
2155         nbv->grp[i].nbl_lists.nnbl = 0;
2156         nbv->grp[i].nbat           = NULL;
2157         nbv->grp[i].kernel_type    = nbnxnkNotSet;
2158
2159         if (i == 0) /* local */
2160         {
2161             pick_nbnxn_kernel(fp, cr, fr->use_simd_kernels,
2162                               nbv->bUseGPU, bEmulateGPU, ir,
2163                               &nbv->grp[i].kernel_type,
2164                               &nbv->grp[i].ewald_excl,
2165                               fr->bNonbonded);
2166         }
2167         else /* non-local */
2168         {
2169             if (nbpu_opt != NULL && strcmp(nbpu_opt, "gpu_cpu") == 0)
2170             {
2171                 /* Use GPU for local, select a CPU kernel for non-local */
2172                 pick_nbnxn_kernel(fp, cr, fr->use_simd_kernels,
2173                                   FALSE, FALSE, ir,
2174                                   &nbv->grp[i].kernel_type,
2175                                   &nbv->grp[i].ewald_excl,
2176                                   fr->bNonbonded);
2177
2178                 bHybridGPURun = TRUE;
2179             }
2180             else
2181             {
2182                 /* Use the same kernel for local and non-local interactions */
2183                 nbv->grp[i].kernel_type = nbv->grp[0].kernel_type;
2184                 nbv->grp[i].ewald_excl  = nbv->grp[0].ewald_excl;
2185             }
2186         }
2187     }
2188
2189     if (nbv->bUseGPU)
2190     {
2191         /* init the NxN GPU data; the last argument tells whether we'll have
2192          * both local and non-local NB calculation on GPU */
2193         nbnxn_cuda_init(fp, &nbv->cu_nbv,
2194                         &fr->hwinfo->gpu_info, fr->gpu_opt,
2195                         cr->rank_pp_intranode,
2196                         (nbv->ngrp > 1) && !bHybridGPURun);
2197
2198         if ((env = getenv("GMX_NB_MIN_CI")) != NULL)
2199         {
2200             char *end;
2201
2202             nbv->min_ci_balanced = strtol(env, &end, 10);
2203             if (!end || (*end != 0) || nbv->min_ci_balanced <= 0)
2204             {
2205                 gmx_fatal(FARGS, "Invalid value passed in GMX_NB_MIN_CI=%s, positive integer required", env);
2206             }
2207
2208             if (debug)
2209             {
2210                 fprintf(debug, "Neighbor-list balancing parameter: %d (passed as env. var.)\n",
2211                         nbv->min_ci_balanced);
2212             }
2213         }
2214         else
2215         {
2216             nbv->min_ci_balanced = nbnxn_cuda_min_ci_balanced(nbv->cu_nbv);
2217             if (debug)
2218             {
2219                 fprintf(debug, "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU multi-processors)\n",
2220                         nbv->min_ci_balanced);
2221             }
2222         }
2223     }
2224     else
2225     {
2226         nbv->min_ci_balanced = 0;
2227     }
2228
2229     *nb_verlet = nbv;
2230
2231     nbnxn_init_search(&nbv->nbs,
2232                       DOMAINDECOMP(cr) ? &cr->dd->nc : NULL,
2233                       DOMAINDECOMP(cr) ? domdec_zones(cr->dd) : NULL,
2234                       bFEP_NonBonded,
2235                       gmx_omp_nthreads_get(emntPairsearch));
2236
2237     for (i = 0; i < nbv->ngrp; i++)
2238     {
2239         if (nbv->grp[0].kernel_type == nbnxnk8x8x8_CUDA)
2240         {
2241             nb_alloc = &pmalloc;
2242             nb_free  = &pfree;
2243         }
2244         else
2245         {
2246             nb_alloc = NULL;
2247             nb_free  = NULL;
2248         }
2249
2250         nbnxn_init_pairlist_set(&nbv->grp[i].nbl_lists,
2251                                 nbnxn_kernel_pairlist_simple(nbv->grp[i].kernel_type),
2252                                 /* 8x8x8 "non-simple" lists are ATM always combined */
2253                                 !nbnxn_kernel_pairlist_simple(nbv->grp[i].kernel_type),
2254                                 nb_alloc, nb_free);
2255
2256         if (i == 0 ||
2257             nbv->grp[0].kernel_type != nbv->grp[i].kernel_type)
2258         {
2259             gmx_bool bSimpleList;
2260             int      enbnxninitcombrule;
2261
2262             bSimpleList = nbnxn_kernel_pairlist_simple(nbv->grp[i].kernel_type);
2263
2264             if (bSimpleList && (fr->vdwtype == evdwCUT && (fr->vdw_modifier == eintmodNONE || fr->vdw_modifier == eintmodPOTSHIFT)))
2265             {
2266                 /* Plain LJ cut-off: we can optimize with combination rules */
2267                 enbnxninitcombrule = enbnxninitcombruleDETECT;
2268             }
2269             else if (fr->vdwtype == evdwPME)
2270             {
2271                 /* LJ-PME: we need to use a combination rule for the grid */
2272                 if (fr->ljpme_combination_rule == eljpmeGEOM)
2273                 {
2274                     enbnxninitcombrule = enbnxninitcombruleGEOM;
2275                 }
2276                 else
2277                 {
2278                     enbnxninitcombrule = enbnxninitcombruleLB;
2279                 }
2280             }
2281             else
2282             {
2283                 /* We use a full combination matrix: no rule required */
2284                 enbnxninitcombrule = enbnxninitcombruleNONE;
2285             }
2286
2287
2288             snew(nbv->grp[i].nbat, 1);
2289             nbnxn_atomdata_init(fp,
2290                                 nbv->grp[i].nbat,
2291                                 nbv->grp[i].kernel_type,
2292                                 enbnxninitcombrule,
2293                                 fr->ntype, fr->nbfp,
2294                                 ir->opts.ngener,
2295                                 bSimpleList ? gmx_omp_nthreads_get(emntNonbonded) : 1,
2296                                 nb_alloc, nb_free);
2297         }
2298         else
2299         {
2300             nbv->grp[i].nbat = nbv->grp[0].nbat;
2301         }
2302     }
2303 }
2304
2305 gmx_bool usingGpu(nonbonded_verlet_t *nbv)
2306 {
2307     return nbv != NULL && nbv->bUseGPU;
2308 }
2309
2310 void init_forcerec(FILE              *fp,
2311                    const output_env_t oenv,
2312                    t_forcerec        *fr,
2313                    t_fcdata          *fcd,
2314                    const t_inputrec  *ir,
2315                    const gmx_mtop_t  *mtop,
2316                    const t_commrec   *cr,
2317                    matrix             box,
2318                    const char        *tabfn,
2319                    const char        *tabafn,
2320                    const char        *tabpfn,
2321                    const char        *tabbfn,
2322                    const char        *nbpu_opt,
2323                    gmx_bool           bNoSolvOpt,
2324                    real               print_force)
2325 {
2326     int            i, m, negp_pp, negptable, egi, egj;
2327     real           rtab;
2328     char          *env;
2329     double         dbl;
2330     const t_block *cgs;
2331     gmx_bool       bGenericKernelOnly;
2332     gmx_bool       bMakeTables, bMakeSeparate14Table, bSomeNormalNbListsAreInUse;
2333     gmx_bool       bFEP_NonBonded;
2334     int           *nm_ind, egp_flags;
2335
2336     if (fr->hwinfo == NULL)
2337     {
2338         /* Detect hardware, gather information.
2339          * In mdrun, hwinfo has already been set before calling init_forcerec.
2340          * Here we ignore GPUs, as tools will not use them anyhow.
2341          */
2342         fr->hwinfo = gmx_detect_hardware(fp, cr, FALSE);
2343     }
2344
2345     /* By default we turn SIMD kernels on, but it might be turned off further down... */
2346     fr->use_simd_kernels = TRUE;
2347
2348     fr->bDomDec = DOMAINDECOMP(cr);
2349
2350     if (check_box(ir->ePBC, box))
2351     {
2352         gmx_fatal(FARGS, check_box(ir->ePBC, box));
2353     }
2354
2355     /* Test particle insertion ? */
2356     if (EI_TPI(ir->eI))
2357     {
2358         /* Set to the size of the molecule to be inserted (the last one) */
2359         /* Because of old style topologies, we have to use the last cg
2360          * instead of the last molecule type.
2361          */
2362         cgs       = &mtop->moltype[mtop->molblock[mtop->nmolblock-1].type].cgs;
2363         fr->n_tpi = cgs->index[cgs->nr] - cgs->index[cgs->nr-1];
2364         if (fr->n_tpi != mtop->mols.index[mtop->mols.nr] - mtop->mols.index[mtop->mols.nr-1])
2365         {
2366             gmx_fatal(FARGS, "The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group.");
2367         }
2368     }
2369     else
2370     {
2371         fr->n_tpi = 0;
2372     }
2373
2374     /* Copy AdResS parameters */
2375     if (ir->bAdress)
2376     {
2377         fr->adress_type           = ir->adress->type;
2378         fr->adress_const_wf       = ir->adress->const_wf;
2379         fr->adress_ex_width       = ir->adress->ex_width;
2380         fr->adress_hy_width       = ir->adress->hy_width;
2381         fr->adress_icor           = ir->adress->icor;
2382         fr->adress_site           = ir->adress->site;
2383         fr->adress_ex_forcecap    = ir->adress->ex_forcecap;
2384         fr->adress_do_hybridpairs = ir->adress->do_hybridpairs;
2385
2386
2387         snew(fr->adress_group_explicit, ir->adress->n_energy_grps);
2388         for (i = 0; i < ir->adress->n_energy_grps; i++)
2389         {
2390             fr->adress_group_explicit[i] = ir->adress->group_explicit[i];
2391         }
2392
2393         fr->n_adress_tf_grps = ir->adress->n_tf_grps;
2394         snew(fr->adress_tf_table_index, fr->n_adress_tf_grps);
2395         for (i = 0; i < fr->n_adress_tf_grps; i++)
2396         {
2397             fr->adress_tf_table_index[i] = ir->adress->tf_table_index[i];
2398         }
2399         copy_rvec(ir->adress->refs, fr->adress_refs);
2400     }
2401     else
2402     {
2403         fr->adress_type           = eAdressOff;
2404         fr->adress_do_hybridpairs = FALSE;
2405     }
2406
2407     /* Copy the user determined parameters */
2408     fr->userint1  = ir->userint1;
2409     fr->userint2  = ir->userint2;
2410     fr->userint3  = ir->userint3;
2411     fr->userint4  = ir->userint4;
2412     fr->userreal1 = ir->userreal1;
2413     fr->userreal2 = ir->userreal2;
2414     fr->userreal3 = ir->userreal3;
2415     fr->userreal4 = ir->userreal4;
2416
2417     /* Shell stuff */
2418     fr->fc_stepsize = ir->fc_stepsize;
2419
2420     /* Free energy */
2421     fr->efep        = ir->efep;
2422     fr->sc_alphavdw = ir->fepvals->sc_alpha;
2423     if (ir->fepvals->bScCoul)
2424     {
2425         fr->sc_alphacoul  = ir->fepvals->sc_alpha;
2426         fr->sc_sigma6_min = pow(ir->fepvals->sc_sigma_min, 6);
2427     }
2428     else
2429     {
2430         fr->sc_alphacoul  = 0;
2431         fr->sc_sigma6_min = 0; /* only needed when bScCoul is on */
2432     }
2433     fr->sc_power      = ir->fepvals->sc_power;
2434     fr->sc_r_power    = ir->fepvals->sc_r_power;
2435     fr->sc_sigma6_def = pow(ir->fepvals->sc_sigma, 6);
2436
2437     env = getenv("GMX_SCSIGMA_MIN");
2438     if (env != NULL)
2439     {
2440         dbl = 0;
2441         sscanf(env, "%20lf", &dbl);
2442         fr->sc_sigma6_min = pow(dbl, 6);
2443         if (fp)
2444         {
2445             fprintf(fp, "Setting the minimum soft core sigma to %g nm\n", dbl);
2446         }
2447     }
2448
2449     fr->bNonbonded = TRUE;
2450     if (getenv("GMX_NO_NONBONDED") != NULL)
2451     {
2452         /* turn off non-bonded calculations */
2453         fr->bNonbonded = FALSE;
2454         md_print_warn(cr, fp,
2455                       "Found environment variable GMX_NO_NONBONDED.\n"
2456                       "Disabling nonbonded calculations.\n");
2457     }
2458
2459     bGenericKernelOnly = FALSE;
2460
2461     /* We now check in the NS code whether a particular combination of interactions
2462      * can be used with water optimization, and disable it if that is not the case.
2463      */
2464
2465     if (getenv("GMX_NB_GENERIC") != NULL)
2466     {
2467         if (fp != NULL)
2468         {
2469             fprintf(fp,
2470                     "Found environment variable GMX_NB_GENERIC.\n"
2471                     "Disabling all interaction-specific nonbonded kernels, will only\n"
2472                     "use the slow generic ones in src/gmxlib/nonbonded/nb_generic.c\n\n");
2473         }
2474         bGenericKernelOnly = TRUE;
2475     }
2476
2477     if (bGenericKernelOnly == TRUE)
2478     {
2479         bNoSolvOpt         = TRUE;
2480     }
2481
2482     if ( (getenv("GMX_DISABLE_SIMD_KERNELS") != NULL) || (getenv("GMX_NOOPTIMIZEDKERNELS") != NULL) )
2483     {
2484         fr->use_simd_kernels = FALSE;
2485         if (fp != NULL)
2486         {
2487             fprintf(fp,
2488                     "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
2489                     "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
2490                     "(e.g. SSE2/SSE4.1/AVX).\n\n");
2491         }
2492     }
2493
2494     fr->bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
2495
2496     /* Check if we can/should do all-vs-all kernels */
2497     fr->bAllvsAll       = can_use_allvsall(ir, FALSE, NULL, NULL);
2498     fr->AllvsAll_work   = NULL;
2499     fr->AllvsAll_workgb = NULL;
2500
2501     /* All-vs-all kernels have not been implemented in 4.6, and
2502      * the SIMD group kernels are also buggy in this case. Non-SIMD
2503      * group kernels are OK. See Redmine #1249. */
2504     if (fr->bAllvsAll)
2505     {
2506         fr->bAllvsAll            = FALSE;
2507         fr->use_simd_kernels     = FALSE;
2508         if (fp != NULL)
2509         {
2510             fprintf(fp,
2511                     "\nYour simulation settings would have triggered the efficient all-vs-all\n"
2512                     "kernels in GROMACS 4.5, but these have not been implemented in GROMACS\n"
2513                     "4.6. Also, we can't use the accelerated SIMD kernels here because\n"
2514                     "of an unfixed bug. The reference C kernels are correct, though, so\n"
2515                     "we are proceeding by disabling all CPU architecture-specific\n"
2516                     "(e.g. SSE2/SSE4/AVX) routines. If performance is important, please\n"
2517                     "use GROMACS 4.5.7 or try cutoff-scheme = Verlet.\n\n");
2518         }
2519     }
2520
2521     /* Neighbour searching stuff */
2522     fr->cutoff_scheme = ir->cutoff_scheme;
2523     fr->bGrid         = (ir->ns_type == ensGRID);
2524     fr->ePBC          = ir->ePBC;
2525
2526     if (fr->cutoff_scheme == ecutsGROUP)
2527     {
2528         const char *note = "NOTE: This file uses the deprecated 'group' cutoff_scheme. This will be\n"
2529             "removed in a future release when 'verlet' supports all interaction forms.\n";
2530
2531         if (MASTER(cr))
2532         {
2533             fprintf(stderr, "\n%s\n", note);
2534         }
2535         if (fp != NULL)
2536         {
2537             fprintf(fp, "\n%s\n", note);
2538         }
2539     }
2540
2541     /* Determine if we will do PBC for distances in bonded interactions */
2542     if (fr->ePBC == epbcNONE)
2543     {
2544         fr->bMolPBC = FALSE;
2545     }
2546     else
2547     {
2548         if (!DOMAINDECOMP(cr))
2549         {
2550             /* The group cut-off scheme and SHAKE assume charge groups
2551              * are whole, but not using molpbc is faster in most cases.
2552              */
2553             if (fr->cutoff_scheme == ecutsGROUP ||
2554                 (ir->eConstrAlg == econtSHAKE &&
2555                  (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 ||
2556                   gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0)))
2557             {
2558                 fr->bMolPBC = ir->bPeriodicMols;
2559             }
2560             else
2561             {
2562                 fr->bMolPBC = TRUE;
2563                 if (getenv("GMX_USE_GRAPH") != NULL)
2564                 {
2565                     fr->bMolPBC = FALSE;
2566                     if (fp)
2567                     {
2568                         fprintf(fp, "\nGMX_MOLPBC is set, using the graph for bonded interactions\n\n");
2569                     }
2570                 }
2571             }
2572         }
2573         else
2574         {
2575             fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->ePBC);
2576         }
2577     }
2578     fr->bGB = (ir->implicit_solvent == eisGBSA);
2579
2580     fr->rc_scaling = ir->refcoord_scaling;
2581     copy_rvec(ir->posres_com, fr->posres_com);
2582     copy_rvec(ir->posres_comB, fr->posres_comB);
2583     fr->rlist                    = cutoff_inf(ir->rlist);
2584     fr->rlistlong                = cutoff_inf(ir->rlistlong);
2585     fr->eeltype                  = ir->coulombtype;
2586     fr->vdwtype                  = ir->vdwtype;
2587     fr->ljpme_combination_rule   = ir->ljpme_combination_rule;
2588
2589     fr->coulomb_modifier = ir->coulomb_modifier;
2590     fr->vdw_modifier     = ir->vdw_modifier;
2591
2592     /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
2593     switch (fr->eeltype)
2594     {
2595         case eelCUT:
2596             fr->nbkernel_elec_interaction = (fr->bGB) ? GMX_NBKERNEL_ELEC_GENERALIZEDBORN : GMX_NBKERNEL_ELEC_COULOMB;
2597             break;
2598
2599         case eelRF:
2600         case eelGRF:
2601         case eelRF_NEC:
2602             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_REACTIONFIELD;
2603             break;
2604
2605         case eelRF_ZERO:
2606             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_REACTIONFIELD;
2607             fr->coulomb_modifier          = eintmodEXACTCUTOFF;
2608             break;
2609
2610         case eelSWITCH:
2611         case eelSHIFT:
2612         case eelUSER:
2613         case eelENCADSHIFT:
2614         case eelPMESWITCH:
2615         case eelPMEUSER:
2616         case eelPMEUSERSWITCH:
2617             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_CUBICSPLINETABLE;
2618             break;
2619
2620         case eelPME:
2621         case eelEWALD:
2622             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_EWALD;
2623             break;
2624
2625         default:
2626             gmx_fatal(FARGS, "Unsupported electrostatic interaction: %s", eel_names[fr->eeltype]);
2627             break;
2628     }
2629
2630     /* Vdw: Translate from mdp settings to kernel format */
2631     switch (fr->vdwtype)
2632     {
2633         case evdwCUT:
2634             if (fr->bBHAM)
2635             {
2636                 fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_BUCKINGHAM;
2637             }
2638             else
2639             {
2640                 fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_LENNARDJONES;
2641             }
2642             break;
2643         case evdwPME:
2644             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_LJEWALD;
2645             break;
2646
2647         case evdwSWITCH:
2648         case evdwSHIFT:
2649         case evdwUSER:
2650         case evdwENCADSHIFT:
2651             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE;
2652             break;
2653
2654         default:
2655             gmx_fatal(FARGS, "Unsupported vdw interaction: %s", evdw_names[fr->vdwtype]);
2656             break;
2657     }
2658
2659     /* These start out identical to ir, but might be altered if we e.g. tabulate the interaction in the kernel */
2660     fr->nbkernel_elec_modifier    = fr->coulomb_modifier;
2661     fr->nbkernel_vdw_modifier     = fr->vdw_modifier;
2662
2663     fr->rvdw             = cutoff_inf(ir->rvdw);
2664     fr->rvdw_switch      = ir->rvdw_switch;
2665     fr->rcoulomb         = cutoff_inf(ir->rcoulomb);
2666     fr->rcoulomb_switch  = ir->rcoulomb_switch;
2667
2668     fr->bTwinRange = fr->rlistlong > fr->rlist;
2669     fr->bEwald     = (EEL_PME(fr->eeltype) || fr->eeltype == eelEWALD);
2670
2671     fr->reppow     = mtop->ffparams.reppow;
2672
2673     if (ir->cutoff_scheme == ecutsGROUP)
2674     {
2675         fr->bvdwtab    = ((fr->vdwtype != evdwCUT || !gmx_within_tol(fr->reppow, 12.0, 10*GMX_DOUBLE_EPS))
2676                           && !EVDW_PME(fr->vdwtype));
2677         /* We have special kernels for standard Ewald and PME, but the pme-switch ones are tabulated above */
2678         fr->bcoultab   = !(fr->eeltype == eelCUT ||
2679                            fr->eeltype == eelEWALD ||
2680                            fr->eeltype == eelPME ||
2681                            fr->eeltype == eelRF ||
2682                            fr->eeltype == eelRF_ZERO);
2683
2684         /* If the user absolutely wants different switch/shift settings for coul/vdw, it is likely
2685          * going to be faster to tabulate the interaction than calling the generic kernel.
2686          * However, if generic kernels have been requested we keep things analytically.
2687          */
2688         if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH &&
2689             fr->nbkernel_vdw_modifier == eintmodPOTSWITCH &&
2690             bGenericKernelOnly == FALSE)
2691         {
2692             if ((fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw))
2693             {
2694                 fr->bcoultab = TRUE;
2695                 /* Once we tabulate electrostatics, we can use the switch function for LJ,
2696                  * which would otherwise need two tables.
2697                  */
2698             }
2699         }
2700         else if ((fr->nbkernel_elec_modifier == eintmodPOTSHIFT && fr->nbkernel_vdw_modifier == eintmodPOTSHIFT) ||
2701                  ((fr->nbkernel_elec_interaction == GMX_NBKERNEL_ELEC_REACTIONFIELD &&
2702                    fr->nbkernel_elec_modifier == eintmodEXACTCUTOFF &&
2703                    (fr->nbkernel_vdw_modifier == eintmodPOTSWITCH || fr->nbkernel_vdw_modifier == eintmodPOTSHIFT))))
2704         {
2705             if ((fr->rcoulomb != fr->rvdw) && (bGenericKernelOnly == FALSE))
2706             {
2707                 fr->bcoultab = TRUE;
2708             }
2709         }
2710
2711         if (fr->nbkernel_elec_modifier == eintmodFORCESWITCH)
2712         {
2713             fr->bcoultab = TRUE;
2714         }
2715         if (fr->nbkernel_vdw_modifier == eintmodFORCESWITCH)
2716         {
2717             fr->bvdwtab = TRUE;
2718         }
2719
2720         if (getenv("GMX_REQUIRE_TABLES"))
2721         {
2722             fr->bvdwtab  = TRUE;
2723             fr->bcoultab = TRUE;
2724         }
2725
2726         if (fp)
2727         {
2728             fprintf(fp, "Table routines are used for coulomb: %s\n", bool_names[fr->bcoultab]);
2729             fprintf(fp, "Table routines are used for vdw:     %s\n", bool_names[fr->bvdwtab ]);
2730         }
2731
2732         if (fr->bvdwtab == TRUE)
2733         {
2734             fr->nbkernel_vdw_interaction = GMX_NBKERNEL_VDW_CUBICSPLINETABLE;
2735             fr->nbkernel_vdw_modifier    = eintmodNONE;
2736         }
2737         if (fr->bcoultab == TRUE)
2738         {
2739             fr->nbkernel_elec_interaction = GMX_NBKERNEL_ELEC_CUBICSPLINETABLE;
2740             fr->nbkernel_elec_modifier    = eintmodNONE;
2741         }
2742     }
2743
2744     if (ir->cutoff_scheme == ecutsVERLET)
2745     {
2746         if (!gmx_within_tol(fr->reppow, 12.0, 10*GMX_DOUBLE_EPS))
2747         {
2748             gmx_fatal(FARGS, "Cut-off scheme %S only supports LJ repulsion power 12", ecutscheme_names[ir->cutoff_scheme]);
2749         }
2750         fr->bvdwtab  = FALSE;
2751         fr->bcoultab = FALSE;
2752     }
2753
2754     /* Tables are used for direct ewald sum */
2755     if (fr->bEwald)
2756     {
2757         if (EEL_PME(ir->coulombtype))
2758         {
2759             if (fp)
2760             {
2761                 fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n");
2762             }
2763             if (ir->coulombtype == eelP3M_AD)
2764             {
2765                 please_cite(fp, "Hockney1988");
2766                 please_cite(fp, "Ballenegger2012");
2767             }
2768             else
2769             {
2770                 please_cite(fp, "Essmann95a");
2771             }
2772
2773             if (ir->ewald_geometry == eewg3DC)
2774             {
2775                 if (fp)
2776                 {
2777                     fprintf(fp, "Using the Ewald3DC correction for systems with a slab geometry.\n");
2778                 }
2779                 please_cite(fp, "In-Chul99a");
2780             }
2781         }
2782         fr->ewaldcoeff_q = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
2783         init_ewald_tab(&(fr->ewald_table), ir, fp);
2784         if (fp)
2785         {
2786             fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n",
2787                     1/fr->ewaldcoeff_q);
2788         }
2789     }
2790
2791     if (EVDW_PME(ir->vdwtype))
2792     {
2793         if (fp)
2794         {
2795             fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n");
2796         }
2797         please_cite(fp, "Essmann95a");
2798         fr->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj);
2799         if (fp)
2800         {
2801             fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n",
2802                     1/fr->ewaldcoeff_lj);
2803         }
2804     }
2805
2806     /* Electrostatics */
2807     fr->epsilon_r       = ir->epsilon_r;
2808     fr->epsilon_rf      = ir->epsilon_rf;
2809     fr->fudgeQQ         = mtop->ffparams.fudgeQQ;
2810
2811     /* Parameters for generalized RF */
2812     fr->zsquare = 0.0;
2813     fr->temp    = 0.0;
2814
2815     if (fr->eeltype == eelGRF)
2816     {
2817         init_generalized_rf(fp, mtop, ir, fr);
2818     }
2819
2820     fr->bF_NoVirSum = (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype) ||
2821                        gmx_mtop_ftype_count(mtop, F_POSRES) > 0 ||
2822                        gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 ||
2823                        IR_ELEC_FIELD(*ir) ||
2824                        (fr->adress_icor != eAdressICOff)
2825                        );
2826
2827     if (fr->cutoff_scheme == ecutsGROUP &&
2828         ncg_mtop(mtop) > fr->cg_nalloc && !DOMAINDECOMP(cr))
2829     {
2830         /* Count the total number of charge groups */
2831         fr->cg_nalloc = ncg_mtop(mtop);
2832         srenew(fr->cg_cm, fr->cg_nalloc);
2833     }
2834     if (fr->shift_vec == NULL)
2835     {
2836         snew(fr->shift_vec, SHIFTS);
2837     }
2838
2839     if (fr->fshift == NULL)
2840     {
2841         snew(fr->fshift, SHIFTS);
2842     }
2843
2844     if (fr->nbfp == NULL)
2845     {
2846         fr->ntype = mtop->ffparams.atnr;
2847         fr->nbfp  = mk_nbfp(&mtop->ffparams, fr->bBHAM);
2848         if (EVDW_PME(fr->vdwtype))
2849         {
2850             fr->ljpme_c6grid  = make_ljpme_c6grid(&mtop->ffparams, fr);
2851         }
2852     }
2853
2854     /* Copy the energy group exclusions */
2855     fr->egp_flags = ir->opts.egp_flags;
2856
2857     /* Van der Waals stuff */
2858     if ((fr->vdwtype != evdwCUT) && (fr->vdwtype != evdwUSER) && !fr->bBHAM)
2859     {
2860         if (fr->rvdw_switch >= fr->rvdw)
2861         {
2862             gmx_fatal(FARGS, "rvdw_switch (%f) must be < rvdw (%f)",
2863                       fr->rvdw_switch, fr->rvdw);
2864         }
2865         if (fp)
2866         {
2867             fprintf(fp, "Using %s Lennard-Jones, switch between %g and %g nm\n",
2868                     (fr->eeltype == eelSWITCH) ? "switched" : "shifted",
2869                     fr->rvdw_switch, fr->rvdw);
2870         }
2871     }
2872
2873     if (fr->bBHAM && EVDW_PME(fr->vdwtype))
2874     {
2875         gmx_fatal(FARGS, "LJ PME not supported with Buckingham");
2876     }
2877
2878     if (fr->bBHAM && (fr->vdwtype == evdwSHIFT || fr->vdwtype == evdwSWITCH))
2879     {
2880         gmx_fatal(FARGS, "Switch/shift interaction not supported with Buckingham");
2881     }
2882
2883     if (fr->bBHAM && fr->cutoff_scheme == ecutsVERLET)
2884     {
2885         gmx_fatal(FARGS, "Verlet cutoff-scheme is not supported with Buckingham");
2886     }
2887
2888     if (fp)
2889     {
2890         fprintf(fp, "Cut-off's:   NS: %g   Coulomb: %g   %s: %g\n",
2891                 fr->rlist, fr->rcoulomb, fr->bBHAM ? "BHAM" : "LJ", fr->rvdw);
2892     }
2893
2894     fr->eDispCorr = ir->eDispCorr;
2895     if (ir->eDispCorr != edispcNO)
2896     {
2897         set_avcsixtwelve(fp, fr, mtop);
2898     }
2899
2900     if (fr->bBHAM)
2901     {
2902         set_bham_b_max(fp, fr, mtop);
2903     }
2904
2905     fr->gb_epsilon_solvent = ir->gb_epsilon_solvent;
2906
2907     /* Copy the GBSA data (radius, volume and surftens for each
2908      * atomtype) from the topology atomtype section to forcerec.
2909      */
2910     snew(fr->atype_radius, fr->ntype);
2911     snew(fr->atype_vol, fr->ntype);
2912     snew(fr->atype_surftens, fr->ntype);
2913     snew(fr->atype_gb_radius, fr->ntype);
2914     snew(fr->atype_S_hct, fr->ntype);
2915
2916     if (mtop->atomtypes.nr > 0)
2917     {
2918         for (i = 0; i < fr->ntype; i++)
2919         {
2920             fr->atype_radius[i] = mtop->atomtypes.radius[i];
2921         }
2922         for (i = 0; i < fr->ntype; i++)
2923         {
2924             fr->atype_vol[i] = mtop->atomtypes.vol[i];
2925         }
2926         for (i = 0; i < fr->ntype; i++)
2927         {
2928             fr->atype_surftens[i] = mtop->atomtypes.surftens[i];
2929         }
2930         for (i = 0; i < fr->ntype; i++)
2931         {
2932             fr->atype_gb_radius[i] = mtop->atomtypes.gb_radius[i];
2933         }
2934         for (i = 0; i < fr->ntype; i++)
2935         {
2936             fr->atype_S_hct[i] = mtop->atomtypes.S_hct[i];
2937         }
2938     }
2939
2940     /* Generate the GB table if needed */
2941     if (fr->bGB)
2942     {
2943 #ifdef GMX_DOUBLE
2944         fr->gbtabscale = 2000;
2945 #else
2946         fr->gbtabscale = 500;
2947 #endif
2948
2949         fr->gbtabr = 100;
2950         fr->gbtab  = make_gb_table(oenv, fr);
2951
2952         init_gb(&fr->born, fr, ir, mtop, ir->gb_algorithm);
2953
2954         /* Copy local gb data (for dd, this is done in dd_partition_system) */
2955         if (!DOMAINDECOMP(cr))
2956         {
2957             make_local_gb(cr, fr->born, ir->gb_algorithm);
2958         }
2959     }
2960
2961     /* Set the charge scaling */
2962     if (fr->epsilon_r != 0)
2963     {
2964         fr->epsfac = ONE_4PI_EPS0/fr->epsilon_r;
2965     }
2966     else
2967     {
2968         /* eps = 0 is infinite dieletric: no coulomb interactions */
2969         fr->epsfac = 0;
2970     }
2971
2972     /* Reaction field constants */
2973     if (EEL_RF(fr->eeltype))
2974     {
2975         calc_rffac(fp, fr->eeltype, fr->epsilon_r, fr->epsilon_rf,
2976                    fr->rcoulomb, fr->temp, fr->zsquare, box,
2977                    &fr->kappa, &fr->k_rf, &fr->c_rf);
2978     }
2979
2980     /*This now calculates sum for q and c6*/
2981     set_chargesum(fp, fr, mtop);
2982
2983     /* if we are using LR electrostatics, and they are tabulated,
2984      * the tables will contain modified coulomb interactions.
2985      * Since we want to use the non-shifted ones for 1-4
2986      * coulombic interactions, we must have an extra set of tables.
2987      */
2988
2989     /* Construct tables.
2990      * A little unnecessary to make both vdw and coul tables sometimes,
2991      * but what the heck... */
2992
2993     bMakeTables = fr->bcoultab || fr->bvdwtab || fr->bEwald ||
2994         (ir->eDispCorr != edispcNO && ir_vdw_switched(ir));
2995
2996     bMakeSeparate14Table = ((!bMakeTables || fr->eeltype != eelCUT || fr->vdwtype != evdwCUT ||
2997                              fr->coulomb_modifier != eintmodNONE ||
2998                              fr->vdw_modifier != eintmodNONE ||
2999                              fr->bBHAM || fr->bEwald) &&
3000                             (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 ||
3001                              gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0 ||
3002                              gmx_mtop_ftype_count(mtop, F_LJC_PAIRS_NB) > 0));
3003
3004     negp_pp   = ir->opts.ngener - ir->nwall;
3005     negptable = 0;
3006     if (!bMakeTables)
3007     {
3008         bSomeNormalNbListsAreInUse = TRUE;
3009         fr->nnblists               = 1;
3010     }
3011     else
3012     {
3013         bSomeNormalNbListsAreInUse = (ir->eDispCorr != edispcNO);
3014         for (egi = 0; egi < negp_pp; egi++)
3015         {
3016             for (egj = egi; egj < negp_pp; egj++)
3017             {
3018                 egp_flags = ir->opts.egp_flags[GID(egi, egj, ir->opts.ngener)];
3019                 if (!(egp_flags & EGP_EXCL))
3020                 {
3021                     if (egp_flags & EGP_TABLE)
3022                     {
3023                         negptable++;
3024                     }
3025                     else
3026                     {
3027                         bSomeNormalNbListsAreInUse = TRUE;
3028                     }
3029                 }
3030             }
3031         }
3032         if (bSomeNormalNbListsAreInUse)
3033         {
3034             fr->nnblists = negptable + 1;
3035         }
3036         else
3037         {
3038             fr->nnblists = negptable;
3039         }
3040         if (fr->nnblists > 1)
3041         {
3042             snew(fr->gid2nblists, ir->opts.ngener*ir->opts.ngener);
3043         }
3044     }
3045
3046     if (ir->adress)
3047     {
3048         fr->nnblists *= 2;
3049     }
3050
3051     snew(fr->nblists, fr->nnblists);
3052
3053     /* This code automatically gives table length tabext without cut-off's,
3054      * in that case grompp should already have checked that we do not need
3055      * normal tables and we only generate tables for 1-4 interactions.
3056      */
3057     rtab = ir->rlistlong + ir->tabext;
3058
3059     if (bMakeTables)
3060     {
3061         /* make tables for ordinary interactions */
3062         if (bSomeNormalNbListsAreInUse)
3063         {
3064             make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn, NULL, NULL, &fr->nblists[0]);
3065             if (ir->adress)
3066             {
3067                 make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn, NULL, NULL, &fr->nblists[fr->nnblists/2]);
3068             }
3069             if (!bMakeSeparate14Table)
3070             {
3071                 fr->tab14 = fr->nblists[0].table_elec_vdw;
3072             }
3073             m = 1;
3074         }
3075         else
3076         {
3077             m = 0;
3078         }
3079         if (negptable > 0)
3080         {
3081             /* Read the special tables for certain energy group pairs */
3082             nm_ind = mtop->groups.grps[egcENER].nm_ind;
3083             for (egi = 0; egi < negp_pp; egi++)
3084             {
3085                 for (egj = egi; egj < negp_pp; egj++)
3086                 {
3087                     egp_flags = ir->opts.egp_flags[GID(egi, egj, ir->opts.ngener)];
3088                     if ((egp_flags & EGP_TABLE) && !(egp_flags & EGP_EXCL))
3089                     {
3090                         if (fr->nnblists > 1)
3091                         {
3092                             fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = m;
3093                         }
3094                         /* Read the table file with the two energy groups names appended */
3095                         make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn,
3096                                         *mtop->groups.grpname[nm_ind[egi]],
3097                                         *mtop->groups.grpname[nm_ind[egj]],
3098                                         &fr->nblists[m]);
3099                         if (ir->adress)
3100                         {
3101                             make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn,
3102                                             *mtop->groups.grpname[nm_ind[egi]],
3103                                             *mtop->groups.grpname[nm_ind[egj]],
3104                                             &fr->nblists[fr->nnblists/2+m]);
3105                         }
3106                         m++;
3107                     }
3108                     else if (fr->nnblists > 1)
3109                     {
3110                         fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = 0;
3111                     }
3112                 }
3113             }
3114         }
3115     }
3116     else if ((fr->eDispCorr != edispcNO) &&
3117              ((fr->vdw_modifier == eintmodPOTSWITCH) ||
3118               (fr->vdw_modifier == eintmodFORCESWITCH) ||
3119               (fr->vdw_modifier == eintmodPOTSHIFT)))
3120     {
3121         /* Tables might not be used for the potential modifier interactions per se, but
3122          * we still need them to evaluate switch/shift dispersion corrections in this case.
3123          */
3124         make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn, NULL, NULL, &fr->nblists[0]);
3125     }
3126
3127     if (bMakeSeparate14Table)
3128     {
3129         /* generate extra tables with plain Coulomb for 1-4 interactions only */
3130         fr->tab14 = make_tables(fp, oenv, fr, MASTER(cr), tabpfn, rtab,
3131                                 GMX_MAKETABLES_14ONLY);
3132     }
3133
3134     /* Read AdResS Thermo Force table if needed */
3135     if (fr->adress_icor == eAdressICThermoForce)
3136     {
3137         /* old todo replace */
3138
3139         if (ir->adress->n_tf_grps > 0)
3140         {
3141             make_adress_tf_tables(fp, oenv, fr, ir, tabfn, mtop, box);
3142
3143         }
3144         else
3145         {
3146             /* load the default table */
3147             snew(fr->atf_tabs, 1);
3148             fr->atf_tabs[DEFAULT_TF_TABLE] = make_atf_table(fp, oenv, fr, tabafn, box);
3149         }
3150     }
3151
3152     /* Wall stuff */
3153     fr->nwall = ir->nwall;
3154     if (ir->nwall && ir->wall_type == ewtTABLE)
3155     {
3156         make_wall_tables(fp, oenv, ir, tabfn, &mtop->groups, fr);
3157     }
3158
3159     if (fcd && tabbfn)
3160     {
3161         fcd->bondtab  = make_bonded_tables(fp,
3162                                            F_TABBONDS, F_TABBONDSNC,
3163                                            mtop, tabbfn, "b");
3164         fcd->angletab = make_bonded_tables(fp,
3165                                            F_TABANGLES, -1,
3166                                            mtop, tabbfn, "a");
3167         fcd->dihtab   = make_bonded_tables(fp,
3168                                            F_TABDIHS, -1,
3169                                            mtop, tabbfn, "d");
3170     }
3171     else
3172     {
3173         if (debug)
3174         {
3175             fprintf(debug, "No fcdata or table file name passed, can not read table, can not do bonded interactions\n");
3176         }
3177     }
3178
3179     /* QM/MM initialization if requested
3180      */
3181     if (ir->bQMMM)
3182     {
3183         fprintf(stderr, "QM/MM calculation requested.\n");
3184     }
3185
3186     fr->bQMMM      = ir->bQMMM;
3187     fr->qr         = mk_QMMMrec();
3188
3189     /* Set all the static charge group info */
3190     fr->cginfo_mb = init_cginfo_mb(fp, mtop, fr, bNoSolvOpt,
3191                                    &bFEP_NonBonded,
3192                                    &fr->bExcl_IntraCGAll_InterCGNone);
3193     if (DOMAINDECOMP(cr))
3194     {
3195         fr->cginfo = NULL;
3196     }
3197     else
3198     {
3199         fr->cginfo = cginfo_expand(mtop->nmolblock, fr->cginfo_mb);
3200     }
3201
3202     if (!DOMAINDECOMP(cr))
3203     {
3204         forcerec_set_ranges(fr, ncg_mtop(mtop), ncg_mtop(mtop),
3205                             mtop->natoms, mtop->natoms, mtop->natoms);
3206     }
3207
3208     fr->print_force = print_force;
3209
3210
3211     /* coarse load balancing vars */
3212     fr->t_fnbf    = 0.;
3213     fr->t_wait    = 0.;
3214     fr->timesteps = 0;
3215
3216     /* Initialize neighbor search */
3217     init_ns(fp, cr, &fr->ns, fr, mtop);
3218
3219     if (cr->duty & DUTY_PP)
3220     {
3221         gmx_nonbonded_setup(fr, bGenericKernelOnly);
3222         /*
3223            if (ir->bAdress)
3224             {
3225                 gmx_setup_adress_kernels(fp,bGenericKernelOnly);
3226             }
3227          */
3228     }
3229
3230     /* Initialize the thread working data for bonded interactions */
3231     init_forcerec_f_threads(fr, mtop->groups.grps[egcENER].nr);
3232
3233     snew(fr->excl_load, fr->nthreads+1);
3234
3235     if (fr->cutoff_scheme == ecutsVERLET)
3236     {
3237         if (ir->rcoulomb != ir->rvdw)
3238         {
3239             gmx_fatal(FARGS, "With Verlet lists rcoulomb and rvdw should be identical");
3240         }
3241
3242         init_nb_verlet(fp, &fr->nbv, bFEP_NonBonded, ir, fr, cr, nbpu_opt);
3243     }
3244
3245     /* fr->ic is used both by verlet and group kernels (to some extent) now */
3246     init_interaction_const(fp, cr, &fr->ic, fr, rtab);
3247
3248     if (ir->eDispCorr != edispcNO)
3249     {
3250         calc_enervirdiff(fp, ir->eDispCorr, fr);
3251     }
3252 }
3253
3254 #define pr_real(fp, r) fprintf(fp, "%s: %e\n",#r, r)
3255 #define pr_int(fp, i)  fprintf((fp), "%s: %d\n",#i, i)
3256 #define pr_bool(fp, b) fprintf((fp), "%s: %s\n",#b, bool_names[b])
3257
3258 void pr_forcerec(FILE *fp, t_forcerec *fr)
3259 {
3260     int i;
3261
3262     pr_real(fp, fr->rlist);
3263     pr_real(fp, fr->rcoulomb);
3264     pr_real(fp, fr->fudgeQQ);
3265     pr_bool(fp, fr->bGrid);
3266     pr_bool(fp, fr->bTwinRange);
3267     /*pr_int(fp,fr->cg0);
3268        pr_int(fp,fr->hcg);*/
3269     for (i = 0; i < fr->nnblists; i++)
3270     {
3271         pr_int(fp, fr->nblists[i].table_elec_vdw.n);
3272     }
3273     pr_real(fp, fr->rcoulomb_switch);
3274     pr_real(fp, fr->rcoulomb);
3275
3276     fflush(fp);
3277 }
3278
3279 void forcerec_set_excl_load(t_forcerec           *fr,
3280                             const gmx_localtop_t *top)
3281 {
3282     const int *ind, *a;
3283     int        t, i, j, ntot, n, ntarget;
3284
3285     ind = top->excls.index;
3286     a   = top->excls.a;
3287
3288     ntot = 0;
3289     for (i = 0; i < top->excls.nr; i++)
3290     {
3291         for (j = ind[i]; j < ind[i+1]; j++)
3292         {
3293             if (a[j] > i)
3294             {
3295                 ntot++;
3296             }
3297         }
3298     }
3299
3300     fr->excl_load[0] = 0;
3301     n                = 0;
3302     i                = 0;
3303     for (t = 1; t <= fr->nthreads; t++)
3304     {
3305         ntarget = (ntot*t)/fr->nthreads;
3306         while (i < top->excls.nr && n < ntarget)
3307         {
3308             for (j = ind[i]; j < ind[i+1]; j++)
3309             {
3310                 if (a[j] > i)
3311                 {
3312                     n++;
3313                 }
3314             }
3315             i++;
3316         }
3317         fr->excl_load[t] = i;
3318     }
3319 }
3320
3321 /* Frees GPU memory and destroys the CUDA context.
3322  *
3323  * Note that this function needs to be called even if GPUs are not used
3324  * in this run because the PME ranks have no knowledge of whether GPUs
3325  * are used or not, but all ranks need to enter the barrier below.
3326  */
3327 void free_gpu_resources(const t_forcerec     *fr,
3328                         const t_commrec      *cr,
3329                         const gmx_gpu_info_t *gpu_info,
3330                         const gmx_gpu_opt_t  *gpu_opt)
3331 {
3332     gmx_bool bIsPPrankUsingGPU;
3333     char     gpu_err_str[STRLEN];
3334
3335     bIsPPrankUsingGPU = (cr->duty & DUTY_PP) && fr && fr->nbv && fr->nbv->bUseGPU;
3336
3337     if (bIsPPrankUsingGPU)
3338     {
3339         /* free nbnxn data in GPU memory */
3340         nbnxn_cuda_free(fr->nbv->cu_nbv);
3341
3342         /* With tMPI we need to wait for all ranks to finish deallocation before
3343          * destroying the context in free_gpu() as some ranks may be sharing
3344          * GPU and context.
3345          * Note: as only PP ranks need to free GPU resources, so it is safe to
3346          * not call the barrier on PME ranks.
3347          */
3348 #ifdef GMX_THREAD_MPI
3349         if (PAR(cr))
3350         {
3351             gmx_barrier(cr);
3352         }
3353 #endif  /* GMX_THREAD_MPI */
3354
3355         /* uninitialize GPU (by destroying the context) */
3356         if (!free_gpu(cr->rank_pp_intranode, gpu_err_str, gpu_info, gpu_opt))
3357         {
3358             gmx_warning("On rank %d failed to free GPU #%d: %s",
3359                         cr->nodeid, get_current_gpu_device_id(), gpu_err_str);
3360         }
3361     }
3362 }