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