src/gromacs/mdlib/forcerec.c

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