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