src/gromacs/mdlib/forcerec.c

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