src/gromacs/mdlib/forcerec.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "macros.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "physics.h"
  48 #include "force.h"
  49 #include "nonbonded.h"
  50 #include "invblock.h"
  51 #include "names.h"
  52 #include "network.h"
  53 #include "pbc.h"
  54 #include "ns.h"
  55 #include "mshift.h"
  56 #include "txtdump.h"
  57 #include "coulomb.h"
  58 #include "mdrun.h"
  59 #include "domdec.h"
  60 #include "partdec.h"
  61 #include "qmmm.h"
  62 #include "copyrite.h"
  63 #include "mtop_util.h"
  64 #include "gmx_detectcpu.h"
  65
  66 #ifdef _MSC_VER
  67 /* MSVC definition for __cpuid() */
  68 #include <intrin.h>
  69 #endif
  70
  71
  72
  73 t_forcerec *mk_forcerec(void)
  74 {
  75   t_forcerec *fr;
  76
  77   snew(fr,1);
  78
  79   return fr;
  80 }
  81
  82 #ifdef DEBUG
  83 static void pr_nbfp(FILE *fp,real *nbfp,gmx_bool bBHAM,int atnr)
  84 {
  85   int i,j;
  86
  87   for(i=0; (i<atnr); i++) {
  88     for(j=0; (j<atnr); j++) {
  89       fprintf(fp,"%2d - %2d",i,j);
  90       if (bBHAM)
  91         fprintf(fp,"  a=%10g, b=%10g, c=%10g\n",BHAMA(nbfp,atnr,i,j),
  92                 BHAMB(nbfp,atnr,i,j),BHAMC(nbfp,atnr,i,j));
  93       else
  94         fprintf(fp,"  c6=%10g, c12=%10g\n",C6(nbfp,atnr,i,j),
  95                 C12(nbfp,atnr,i,j));
  96     }
  97   }
  98 }
  99 #endif
 100
 101 static real *mk_nbfp(const gmx_ffparams_t *idef,gmx_bool bBHAM)
 102 {
 103   real *nbfp;
 104   int  i,j,k,atnr;
 105
 106   atnr=idef->atnr;
 107   if (bBHAM) {
 108     snew(nbfp,3*atnr*atnr);
 109     for(i=k=0; (i<atnr); i++) {
 110       for(j=0; (j<atnr); j++,k++) {
 111         BHAMA(nbfp,atnr,i,j) = idef->iparams[k].bham.a;
 112         BHAMB(nbfp,atnr,i,j) = idef->iparams[k].bham.b;
 113         BHAMC(nbfp,atnr,i,j) = idef->iparams[k].bham.c;
 114       }
 115     }
 116   }
 117   else {
 118     snew(nbfp,2*atnr*atnr);
 119     for(i=k=0; (i<atnr); i++) {
 120       for(j=0; (j<atnr); j++,k++) {
 121         C6(nbfp,atnr,i,j)   = idef->iparams[k].lj.c6;
 122         C12(nbfp,atnr,i,j)  = idef->iparams[k].lj.c12;
 123       }
 124     }
 125   }
 126   return nbfp;
 127 }
 128
 129 /* This routine sets fr->solvent_opt to the most common solvent in the
 130  * system, e.g. esolSPC or esolTIP4P. It will also mark each charge group in
 131  * the fr->solvent_type array with the correct type (or esolNO).
 132  *
 133  * Charge groups that fulfill the conditions but are not identical to the
 134  * most common one will be marked as esolNO in the solvent_type array.
 135  *
 136  * TIP3p is identical to SPC for these purposes, so we call it
 137  * SPC in the arrays (Apologies to Bill Jorgensen ;-)
 138  *
 139  * NOTE: QM particle should not
 140  * become an optimized solvent. Not even if there is only one charge
 141  * group in the Qm
 142  */
 143
 144 typedef struct
 145 {
 146     int    model;
 147     int    count;
 148     int    vdwtype[4];
 149     real   charge[4];
 150 } solvent_parameters_t;
 151
 152 static void
 153 check_solvent_cg(const gmx_moltype_t   *molt,
 154                  int                   cg0,
 155                  int                   nmol,
 156                  const unsigned char   *qm_grpnr,
 157                  const t_grps          *qm_grps,
 158                  t_forcerec *          fr,
 159                  int                   *n_solvent_parameters,
 160                  solvent_parameters_t  **solvent_parameters_p,
 161                  int                   cginfo,
 162                  int                   *cg_sp)
 163 {
 164     const t_blocka *  excl;
 165     t_atom            *atom;
 166     int               j,k;
 167     int               j0,j1,nj;
 168     gmx_bool              perturbed;
 169     gmx_bool              has_vdw[4];
 170     gmx_bool              match;
 171     real              tmp_charge[4];
 172     int               tmp_vdwtype[4];
 173     int               tjA;
 174     gmx_bool              qm;
 175     solvent_parameters_t *solvent_parameters;
 176
 177     /* We use a list with parameters for each solvent type.
 178      * Every time we discover a new molecule that fulfills the basic
 179      * conditions for a solvent we compare with the previous entries
 180      * in these lists. If the parameters are the same we just increment
 181      * the counter for that type, and otherwise we create a new type
 182      * based on the current molecule.
 183      *
 184      * Once we've finished going through all molecules we check which
 185      * solvent is most common, and mark all those molecules while we
 186      * clear the flag on all others.
 187      */
 188
 189     solvent_parameters = *solvent_parameters_p;
 190
 191     /* Mark the cg first as non optimized */
 192     *cg_sp = -1;
 193
 194     /* Check if this cg has no exclusions with atoms in other charge groups
 195      * and all atoms inside the charge group excluded.
 196      * We only have 3 or 4 atom solvent loops.
 197      */
 198     if (GET_CGINFO_EXCL_INTER(cginfo) ||
 199         !GET_CGINFO_EXCL_INTRA(cginfo))
 200     {
 201         return;
 202     }
 203
 204     /* Get the indices of the first atom in this charge group */
 205     j0     = molt->cgs.index[cg0];
 206     j1     = molt->cgs.index[cg0+1];
 207
 208     /* Number of atoms in our molecule */
 209     nj     = j1 - j0;
 210
 211     if (debug) {
 212         fprintf(debug,
 213                 "Moltype '%s': there are %d atoms in this charge group\n",
 214                 *molt->name,nj);
 215     }
 216
 217     /* Check if it could be an SPC (3 atoms) or TIP4p (4) water,
 218      * otherwise skip it.
 219      */
 220     if (nj<3 || nj>4)
 221     {
 222         return;
 223     }
 224
 225     /* Check if we are doing QM on this group */
 226     qm = FALSE;
 227     if (qm_grpnr != NULL)
 228     {
 229         for(j=j0 ; j<j1 && !qm; j++)
 230         {
 231             qm = (qm_grpnr[j] < qm_grps->nr - 1);
 232         }
 233     }
 234     /* Cannot use solvent optimization with QM */
 235     if (qm)
 236     {
 237         return;
 238     }
 239
 240     atom = molt->atoms.atom;
 241
 242     /* Still looks like a solvent, time to check parameters */
 243
 244     /* If it is perturbed (free energy) we can't use the solvent loops,
 245      * so then we just skip to the next molecule.
 246      */
 247     perturbed = FALSE;
 248
 249     for(j=j0; j<j1 && !perturbed; j++)
 250     {
 251         perturbed = PERTURBED(atom[j]);
 252     }
 253
 254     if (perturbed)
 255     {
 256         return;
 257     }
 258
 259     /* Now it's only a question if the VdW and charge parameters
 260      * are OK. Before doing the check we compare and see if they are
 261      * identical to a possible previous solvent type.
 262      * First we assign the current types and charges.
 263      */
 264     for(j=0; j<nj; j++)
 265     {
 266         tmp_vdwtype[j] = atom[j0+j].type;
 267         tmp_charge[j]  = atom[j0+j].q;
 268     }
 269
 270     /* Does it match any previous solvent type? */
 271     for(k=0 ; k<*n_solvent_parameters; k++)
 272     {
 273         match = TRUE;
 274
 275
 276         /* We can only match SPC with 3 atoms and TIP4p with 4 atoms */
 277         if( (solvent_parameters[k].model==esolSPC   && nj!=3)  ||
 278             (solvent_parameters[k].model==esolTIP4P && nj!=4) )
 279             match = FALSE;
 280
 281         /* Check that types & charges match for all atoms in molecule */
 282         for(j=0 ; j<nj && match==TRUE; j++)
 283         {
 284             if (tmp_vdwtype[j] != solvent_parameters[k].vdwtype[j])
 285             {
 286                 match = FALSE;
 287             }
 288             if(tmp_charge[j] != solvent_parameters[k].charge[j])
 289             {
 290                 match = FALSE;
 291             }
 292         }
 293         if (match == TRUE)
 294         {
 295             /* Congratulations! We have a matched solvent.
 296              * Flag it with this type for later processing.
 297              */
 298             *cg_sp = k;
 299             solvent_parameters[k].count += nmol;
 300
 301             /* We are done with this charge group */
 302             return;
 303         }
 304     }
 305
 306     /* If we get here, we have a tentative new solvent type.
 307      * Before we add it we must check that it fulfills the requirements
 308      * of the solvent optimized loops. First determine which atoms have
 309      * VdW interactions.
 310      */
 311     for(j=0; j<nj; j++)
 312     {
 313         has_vdw[j] = FALSE;
 314         tjA        = tmp_vdwtype[j];
 315
 316         /* Go through all other tpes and see if any have non-zero
 317          * VdW parameters when combined with this one.
 318          */
 319         for(k=0; k<fr->ntype && (has_vdw[j]==FALSE); k++)
 320         {
 321             /* We already checked that the atoms weren't perturbed,
 322              * so we only need to check state A now.
 323              */
 324             if (fr->bBHAM)
 325             {
 326                 has_vdw[j] = (has_vdw[j] ||
 327                               (BHAMA(fr->nbfp,fr->ntype,tjA,k) != 0.0) ||
 328                               (BHAMB(fr->nbfp,fr->ntype,tjA,k) != 0.0) ||
 329                               (BHAMC(fr->nbfp,fr->ntype,tjA,k) != 0.0));
 330             }
 331             else
 332             {
 333                 /* Standard LJ */
 334                 has_vdw[j] = (has_vdw[j] ||
 335                               (C6(fr->nbfp,fr->ntype,tjA,k)  != 0.0) ||
 336                               (C12(fr->nbfp,fr->ntype,tjA,k) != 0.0));
 337             }
 338         }
 339     }
 340
 341     /* Now we know all we need to make the final check and assignment. */
 342     if (nj == 3)
 343     {
 344         /* So, is it an SPC?
 345          * For this we require thatn all atoms have charge,
 346          * the charges on atom 2 & 3 should be the same, and only
 347          * atom 1 should have VdW.
 348          */
 349         if (has_vdw[0] == TRUE &&
 350             has_vdw[1] == FALSE &&
 351             has_vdw[2] == FALSE &&
 352             tmp_charge[0]  != 0 &&
 353             tmp_charge[1]  != 0 &&
 354             tmp_charge[2]  == tmp_charge[1])
 355         {
 356             srenew(solvent_parameters,*n_solvent_parameters+1);
 357             solvent_parameters[*n_solvent_parameters].model = esolSPC;
 358             solvent_parameters[*n_solvent_parameters].count = nmol;
 359             for(k=0;k<3;k++)
 360             {
 361                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 362                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 363             }
 364
 365             *cg_sp = *n_solvent_parameters;
 366             (*n_solvent_parameters)++;
 367         }
 368     }
 369     else if (nj==4)
 370     {
 371         /* Or could it be a TIP4P?
 372          * For this we require thatn atoms 2,3,4 have charge, but not atom 1.
 373          * Only atom 1 should have VdW.
 374          */
 375         if(has_vdw[0] == TRUE &&
 376            has_vdw[1] == FALSE &&
 377            has_vdw[2] == FALSE &&
 378            has_vdw[3] == FALSE &&
 379            tmp_charge[0]  == 0 &&
 380            tmp_charge[1]  != 0 &&
 381            tmp_charge[2]  == tmp_charge[1] &&
 382            tmp_charge[3]  != 0)
 383         {
 384             srenew(solvent_parameters,*n_solvent_parameters+1);
 385             solvent_parameters[*n_solvent_parameters].model = esolTIP4P;
 386             solvent_parameters[*n_solvent_parameters].count = nmol;
 387             for(k=0;k<4;k++)
 388             {
 389                 solvent_parameters[*n_solvent_parameters].vdwtype[k] = tmp_vdwtype[k];
 390                 solvent_parameters[*n_solvent_parameters].charge[k]  = tmp_charge[k];
 391             }
 392
 393             *cg_sp = *n_solvent_parameters;
 394             (*n_solvent_parameters)++;
 395         }
 396     }
 397
 398     *solvent_parameters_p = solvent_parameters;
 399 }
 400
 401 static void
 402 check_solvent(FILE *                fp,
 403               const gmx_mtop_t *    mtop,
 404               t_forcerec *          fr,
 405               cginfo_mb_t           *cginfo_mb)
 406 {
 407     const t_block *   cgs;
 408     const t_block *   mols;
 409     const gmx_moltype_t *molt;
 410     int               mb,mol,cg_mol,at_offset,cg_offset,am,cgm,i,nmol_ch,nmol;
 411     int               n_solvent_parameters;
 412     solvent_parameters_t *solvent_parameters;
 413     int               **cg_sp;
 414     int               bestsp,bestsol;
 415
 416     if (debug)
 417     {
 418         fprintf(debug,"Going to determine what solvent types we have.\n");
 419     }
 420
 421     mols = &mtop->mols;
 422
 423     n_solvent_parameters = 0;
 424     solvent_parameters = NULL;
 425     /* Allocate temporary array for solvent type */
 426     snew(cg_sp,mtop->nmolblock);
 427
 428     cg_offset = 0;
 429     at_offset = 0;
 430     for(mb=0; mb<mtop->nmolblock; mb++)
 431     {
 432         molt = &mtop->moltype[mtop->molblock[mb].type];
 433         cgs  = &molt->cgs;
 434         /* Here we have to loop over all individual molecules
 435          * because we need to check for QMMM particles.
 436          */
 437         snew(cg_sp[mb],cginfo_mb[mb].cg_mod);
 438         nmol_ch = cginfo_mb[mb].cg_mod/cgs->nr;
 439         nmol    = mtop->molblock[mb].nmol/nmol_ch;
 440         for(mol=0; mol<nmol_ch; mol++)
 441         {
 442             cgm = mol*cgs->nr;
 443             am  = mol*cgs->index[cgs->nr];
 444             for(cg_mol=0; cg_mol<cgs->nr; cg_mol++)
 445             {
 446                 check_solvent_cg(molt,cg_mol,nmol,
 447                                  mtop->groups.grpnr[egcQMMM] ?
 448                                  mtop->groups.grpnr[egcQMMM]+at_offset+am : 0,
 449                                  &mtop->groups.grps[egcQMMM],
 450                                  fr,
 451                                  &n_solvent_parameters,&solvent_parameters,
 452                                  cginfo_mb[mb].cginfo[cgm+cg_mol],
 453                                  &cg_sp[mb][cgm+cg_mol]);
 454             }
 455         }
 456         cg_offset += cgs->nr;
 457         at_offset += cgs->index[cgs->nr];
 458     }
 459
 460     /* Puh! We finished going through all charge groups.
 461      * Now find the most common solvent model.
 462      */
 463
 464     /* Most common solvent this far */
 465     bestsp = -2;
 466     for(i=0;i<n_solvent_parameters;i++)
 467     {
 468         if (bestsp == -2 ||
 469             solvent_parameters[i].count > solvent_parameters[bestsp].count)
 470         {
 471             bestsp = i;
 472         }
 473     }
 474
 475     if (bestsp >= 0)
 476     {
 477         bestsol = solvent_parameters[bestsp].model;
 478     }
 479     else
 480     {
 481         bestsol = esolNO;
 482     }
 483
 484 #ifdef DISABLE_WATER_NLIST
 485         bestsol = esolNO;
 486 #endif
 487
 488     fr->nWatMol = 0;
 489     for(mb=0; mb<mtop->nmolblock; mb++)
 490     {
 491         cgs = &mtop->moltype[mtop->molblock[mb].type].cgs;
 492         nmol = (mtop->molblock[mb].nmol*cgs->nr)/cginfo_mb[mb].cg_mod;
 493         for(i=0; i<cginfo_mb[mb].cg_mod; i++)
 494         {
 495             if (cg_sp[mb][i] == bestsp)
 496             {
 497                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i],bestsol);
 498                 fr->nWatMol += nmol;
 499             }
 500             else
 501             {
 502                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[i],esolNO);
 503             }
 504         }
 505         sfree(cg_sp[mb]);
 506     }
 507     sfree(cg_sp);
 508
 509     if (bestsol != esolNO && fp!=NULL)
 510     {
 511         fprintf(fp,"\nEnabling %s-like water optimization for %d molecules.\n\n",
 512                 esol_names[bestsol],
 513                 solvent_parameters[bestsp].count);
 514     }
 515
 516     sfree(solvent_parameters);
 517     fr->solvent_opt = bestsol;
 518 }
 519
 520 static cginfo_mb_t *init_cginfo_mb(FILE *fplog,const gmx_mtop_t *mtop,
 521                                    t_forcerec *fr,gmx_bool bNoSolvOpt)
 522 {
 523     const t_block *cgs;
 524     const t_blocka *excl;
 525     const gmx_moltype_t *molt;
 526     const gmx_molblock_t *molb;
 527     cginfo_mb_t *cginfo_mb;
 528     int  *cginfo;
 529     int  cg_offset,a_offset,cgm,am;
 530     int  mb,m,ncg_tot,cg,a0,a1,gid,ai,j,aj,excl_nalloc;
 531     gmx_bool bId,*bExcl,bExclIntraAll,bExclInter;
 532
 533     ncg_tot = ncg_mtop(mtop);
 534     snew(cginfo_mb,mtop->nmolblock);
 535
 536     excl_nalloc = 10;
 537     snew(bExcl,excl_nalloc);
 538     cg_offset = 0;
 539     a_offset  = 0;
 540     for(mb=0; mb<mtop->nmolblock; mb++)
 541     {
 542         molb = &mtop->molblock[mb];
 543         molt = &mtop->moltype[molb->type];
 544         cgs  = &molt->cgs;
 545         excl = &molt->excls;
 546
 547         /* Check if the cginfo is identical for all molecules in this block.
 548          * If so, we only need an array of the size of one molecule.
 549          * Otherwise we make an array of #mol times #cgs per molecule.
 550          */
 551         bId = TRUE;
 552         am = 0;
 553         for(m=0; m<molb->nmol; m++)
 554         {
 555             am = m*cgs->index[cgs->nr];
 556             for(cg=0; cg<cgs->nr; cg++)
 557             {
 558                 a0 = cgs->index[cg];
 559                 a1 = cgs->index[cg+1];
 560                 if (ggrpnr(&mtop->groups,egcENER,a_offset+am+a0) !=
 561                     ggrpnr(&mtop->groups,egcENER,a_offset   +a0))
 562                 {
 563                     bId = FALSE;
 564                 }
 565                 if (mtop->groups.grpnr[egcQMMM] != NULL)
 566                 {
 567                     for(ai=a0; ai<a1; ai++)
 568                     {
 569                         if (mtop->groups.grpnr[egcQMMM][a_offset+am+ai] !=
 570                             mtop->groups.grpnr[egcQMMM][a_offset   +ai])
 571                         {
 572                             bId = FALSE;
 573                         }
 574                     }
 575                 }
 576             }
 577         }
 578
 579         cginfo_mb[mb].cg_start = cg_offset;
 580         cginfo_mb[mb].cg_end   = cg_offset + molb->nmol*cgs->nr;
 581         cginfo_mb[mb].cg_mod   = (bId ? 1 : molb->nmol)*cgs->nr;
 582         snew(cginfo_mb[mb].cginfo,cginfo_mb[mb].cg_mod);
 583         cginfo = cginfo_mb[mb].cginfo;
 584
 585         for(m=0; m<(bId ? 1 : molb->nmol); m++)
 586         {
 587             cgm = m*cgs->nr;
 588             am  = m*cgs->index[cgs->nr];
 589             for(cg=0; cg<cgs->nr; cg++)
 590             {
 591                 a0 = cgs->index[cg];
 592                 a1 = cgs->index[cg+1];
 593
 594                 /* Store the energy group in cginfo */
 595                 gid = ggrpnr(&mtop->groups,egcENER,a_offset+am+a0);
 596                 SET_CGINFO_GID(cginfo[cgm+cg],gid);
 597
 598                 /* Check the intra/inter charge group exclusions */
 599                 if (a1-a0 > excl_nalloc) {
 600                     excl_nalloc = a1 - a0;
 601                     srenew(bExcl,excl_nalloc);
 602                 }
 603                 /* bExclIntraAll: all intra cg interactions excluded
 604                  * bExclInter:    any inter cg interactions excluded
 605                  */
 606                 bExclIntraAll = TRUE;
 607                 bExclInter    = FALSE;
 608                 for(ai=a0; ai<a1; ai++) {
 609                     /* Clear the exclusion list for atom ai */
 610                     for(aj=a0; aj<a1; aj++) {
 611                         bExcl[aj-a0] = FALSE;
 612                     }
 613                     /* Loop over all the exclusions of atom ai */
 614                     for(j=excl->index[ai]; j<excl->index[ai+1]; j++)
 615                     {
 616                         aj = excl->a[j];
 617                         if (aj < a0 || aj >= a1)
 618                         {
 619                             bExclInter = TRUE;
 620                         }
 621                         else
 622                         {
 623                             bExcl[aj-a0] = TRUE;
 624                         }
 625                     }
 626                     /* Check if ai excludes a0 to a1 */
 627                     for(aj=a0; aj<a1; aj++)
 628                     {
 629                         if (!bExcl[aj-a0])
 630                         {
 631                             bExclIntraAll = FALSE;
 632                         }
 633                     }
 634                 }
 635                 if (bExclIntraAll)
 636                 {
 637                     SET_CGINFO_EXCL_INTRA(cginfo[cgm+cg]);
 638                 }
 639                 if (bExclInter)
 640                 {
 641                     SET_CGINFO_EXCL_INTER(cginfo[cgm+cg]);
 642                 }
 643                 if (a1 - a0 > MAX_CHARGEGROUP_SIZE)
 644                 {
 645                     /* The size in cginfo is currently only read with DD */
 646                     gmx_fatal(FARGS,"A charge group has size %d which is larger than the limit of %d atoms",a1-a0,MAX_CHARGEGROUP_SIZE);
 647                 }
 648                 SET_CGINFO_NATOMS(cginfo[cgm+cg],a1-a0);
 649             }
 650         }
 651         cg_offset += molb->nmol*cgs->nr;
 652         a_offset  += molb->nmol*cgs->index[cgs->nr];
 653     }
 654     sfree(bExcl);
 655
 656     /* the solvent optimizer is called after the QM is initialized,
 657      * because we don't want to have the QM subsystemto become an
 658      * optimized solvent
 659      */
 660
 661     check_solvent(fplog,mtop,fr,cginfo_mb);
 662
 663     if (getenv("GMX_NO_SOLV_OPT"))
 664     {
 665         if (fplog)
 666         {
 667             fprintf(fplog,"Found environment variable GMX_NO_SOLV_OPT.\n"
 668                     "Disabling all solvent optimization\n");
 669         }
 670         fr->solvent_opt = esolNO;
 671     }
 672     if (bNoSolvOpt)
 673     {
 674         fr->solvent_opt = esolNO;
 675     }
 676     if (!fr->solvent_opt)
 677     {
 678         for(mb=0; mb<mtop->nmolblock; mb++)
 679         {
 680             for(cg=0; cg<cginfo_mb[mb].cg_mod; cg++)
 681             {
 682                 SET_CGINFO_SOLOPT(cginfo_mb[mb].cginfo[cg],esolNO);
 683             }
 684         }
 685     }
 686
 687     return cginfo_mb;
 688 }
 689
 690 static int *cginfo_expand(int nmb,cginfo_mb_t *cgi_mb)
 691 {
 692     int ncg,mb,cg;
 693     int *cginfo;
 694
 695     ncg = cgi_mb[nmb-1].cg_end;
 696     snew(cginfo,ncg);
 697     mb = 0;
 698     for(cg=0; cg<ncg; cg++)
 699     {
 700         while (cg >= cgi_mb[mb].cg_end)
 701         {
 702             mb++;
 703         }
 704         cginfo[cg] =
 705             cgi_mb[mb].cginfo[(cg - cgi_mb[mb].cg_start) % cgi_mb[mb].cg_mod];
 706     }
 707
 708     return cginfo;
 709 }
 710
 711 static void set_chargesum(FILE *log,t_forcerec *fr,const gmx_mtop_t *mtop)
 712 {
 713     double qsum;
 714     int    mb,nmol,i;
 715     const t_atoms *atoms;
 716
 717     qsum = 0;
 718     for(mb=0; mb<mtop->nmolblock; mb++)
 719     {
 720         nmol  = mtop->molblock[mb].nmol;
 721         atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 722         for(i=0; i<atoms->nr; i++)
 723         {
 724             qsum += nmol*atoms->atom[i].q;
 725         }
 726     }
 727     fr->qsum[0] = qsum;
 728     if (fr->efep != efepNO)
 729     {
 730         qsum = 0;
 731         for(mb=0; mb<mtop->nmolblock; mb++)
 732         {
 733             nmol  = mtop->molblock[mb].nmol;
 734             atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 735             for(i=0; i<atoms->nr; i++)
 736             {
 737                 qsum += nmol*atoms->atom[i].qB;
 738             }
 739             fr->qsum[1] = qsum;
 740         }
 741     }
 742     else
 743     {
 744         fr->qsum[1] = fr->qsum[0];
 745     }
 746     if (log) {
 747         if (fr->efep == efepNO)
 748             fprintf(log,"System total charge: %.3f\n",fr->qsum[0]);
 749         else
 750             fprintf(log,"System total charge, top. A: %.3f top. B: %.3f\n",
 751                     fr->qsum[0],fr->qsum[1]);
 752     }
 753 }
 754
 755 void update_forcerec(FILE *log,t_forcerec *fr,matrix box)
 756 {
 757     if (fr->eeltype == eelGRF)
 758     {
 759         calc_rffac(NULL,fr->eeltype,fr->epsilon_r,fr->epsilon_rf,
 760                    fr->rcoulomb,fr->temp,fr->zsquare,box,
 761                    &fr->kappa,&fr->k_rf,&fr->c_rf);
 762     }
 763 }
 764
 765 void set_avcsixtwelve(FILE *fplog,t_forcerec *fr,const gmx_mtop_t *mtop)
 766 {
 767     const t_atoms *atoms,*atoms_tpi;
 768     const t_blocka *excl;
 769     int    mb,nmol,nmolc,i,j,tpi,tpj,j1,j2,k,n,nexcl,q;
 770 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
 771     long long int  npair,npair_ij,tmpi,tmpj;
 772 #else
 773     double npair, npair_ij,tmpi,tmpj;
 774 #endif
 775     double csix,ctwelve;
 776     int    ntp,*typecount;
 777     gmx_bool   bBHAM;
 778     real   *nbfp;
 779
 780     ntp = fr->ntype;
 781     bBHAM = fr->bBHAM;
 782     nbfp = fr->nbfp;
 783
 784     for(q=0; q<(fr->efep==efepNO ? 1 : 2); q++) {
 785         csix = 0;
 786         ctwelve = 0;
 787         npair = 0;
 788         nexcl = 0;
 789         if (!fr->n_tpi) {
 790             /* Count the types so we avoid natoms^2 operations */
 791             snew(typecount,ntp);
 792             for(mb=0; mb<mtop->nmolblock; mb++) {
 793                 nmol  = mtop->molblock[mb].nmol;
 794                 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 795                 for(i=0; i<atoms->nr; i++) {
 796                     if (q == 0)
 797                     {
 798                         tpi = atoms->atom[i].type;
 799                     }
 800                     else
 801                     {
 802                         tpi = atoms->atom[i].typeB;
 803                     }
 804                     typecount[tpi] += nmol;
 805                 }
 806             }
 807             for(tpi=0; tpi<ntp; tpi++) {
 808                 for(tpj=tpi; tpj<ntp; tpj++) {
 809                     tmpi = typecount[tpi];
 810                     tmpj = typecount[tpj];
 811                     if (tpi != tpj)
 812                     {
 813                         npair_ij = tmpi*tmpj;
 814                     }
 815                     else
 816                     {
 817                         npair_ij = tmpi*(tmpi - 1)/2;
 818                     }
 819                     if (bBHAM) {
 820                         csix    += npair_ij*BHAMC(nbfp,ntp,tpi,tpj);
 821                     } else {
 822                         csix    += npair_ij*   C6(nbfp,ntp,tpi,tpj);
 823                         ctwelve += npair_ij*  C12(nbfp,ntp,tpi,tpj);
 824                     }
 825                     npair += npair_ij;
 826                 }
 827             }
 828             sfree(typecount);
 829             /* Subtract the excluded pairs.
 830              * The main reason for substracting exclusions is that in some cases
 831              * some combinations might never occur and the parameters could have
 832              * any value. These unused values should not influence the dispersion
 833              * correction.
 834              */
 835             for(mb=0; mb<mtop->nmolblock; mb++) {
 836                 nmol  = mtop->molblock[mb].nmol;
 837                 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 838                 excl  = &mtop->moltype[mtop->molblock[mb].type].excls;
 839                 for(i=0; (i<atoms->nr); i++) {
 840                     if (q == 0)
 841                     {
 842                         tpi = atoms->atom[i].type;
 843                     }
 844                     else
 845                     {
 846                         tpi = atoms->atom[i].typeB;
 847                     }
 848                     j1  = excl->index[i];
 849                     j2  = excl->index[i+1];
 850                     for(j=j1; j<j2; j++) {
 851                         k = excl->a[j];
 852                         if (k > i)
 853                         {
 854                             if (q == 0)
 855                             {
 856                                 tpj = atoms->atom[k].type;
 857                             }
 858                             else
 859                             {
 860                                 tpj = atoms->atom[k].typeB;
 861                             }
 862                             if (bBHAM) {
 863                                csix -= nmol*BHAMC(nbfp,ntp,tpi,tpj);
 864                             } else {
 865                                 csix    -= nmol*C6 (nbfp,ntp,tpi,tpj);
 866                                 ctwelve -= nmol*C12(nbfp,ntp,tpi,tpj);
 867                             }
 868                             nexcl += nmol;
 869                         }
 870                     }
 871                 }
 872             }
 873         } else {
 874             /* Only correct for the interaction of the test particle
 875              * with the rest of the system.
 876              */
 877             atoms_tpi =
 878                 &mtop->moltype[mtop->molblock[mtop->nmolblock-1].type].atoms;
 879
 880             npair = 0;
 881             for(mb=0; mb<mtop->nmolblock; mb++) {
 882                 nmol  = mtop->molblock[mb].nmol;
 883                 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 884                 for(j=0; j<atoms->nr; j++) {
 885                     nmolc = nmol;
 886                     /* Remove the interaction of the test charge group
 887                      * with itself.
 888                      */
 889                     if (mb == mtop->nmolblock-1)
 890                     {
 891                         nmolc--;
 892
 893                         if (mb == 0 && nmol == 1)
 894                         {
 895                             gmx_fatal(FARGS,"Old format tpr with TPI, please generate a new tpr file");
 896                         }
 897                     }
 898                     if (q == 0)
 899                     {
 900                         tpj = atoms->atom[j].type;
 901                     }
 902                     else
 903                     {
 904                         tpj = atoms->atom[j].typeB;
 905                     }
 906                     for(i=0; i<fr->n_tpi; i++)
 907                     {
 908                         if (q == 0)
 909                         {
 910                             tpi = atoms_tpi->atom[i].type;
 911                         }
 912                         else
 913                         {
 914                             tpi = atoms_tpi->atom[i].typeB;
 915                         }
 916                         if (bBHAM)
 917                         {
 918                             csix    += nmolc*BHAMC(nbfp,ntp,tpi,tpj);
 919                         }
 920                         else
 921                         {
 922                             csix    += nmolc*C6 (nbfp,ntp,tpi,tpj);
 923                             ctwelve += nmolc*C12(nbfp,ntp,tpi,tpj);
 924                         }
 925                         npair += nmolc;
 926                     }
 927                 }
 928             }
 929         }
 930         if (npair - nexcl <= 0 && fplog) {
 931             fprintf(fplog,"\nWARNING: There are no atom pairs for dispersion correction\n\n");
 932             csix     = 0;
 933             ctwelve  = 0;
 934         } else {
 935             csix    /= npair - nexcl;
 936             ctwelve /= npair - nexcl;
 937         }
 938         if (debug) {
 939             fprintf(debug,"Counted %d exclusions\n",nexcl);
 940             fprintf(debug,"Average C6 parameter is: %10g\n",(double)csix);
 941             fprintf(debug,"Average C12 parameter is: %10g\n",(double)ctwelve);
 942         }
 943         fr->avcsix[q]    = csix;
 944         fr->avctwelve[q] = ctwelve;
 945     }
 946     if (fplog != NULL)
 947     {
 948         if (fr->eDispCorr == edispcAllEner ||
 949             fr->eDispCorr == edispcAllEnerPres)
 950         {
 951             fprintf(fplog,"Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
 952                     fr->avcsix[0],fr->avctwelve[0]);
 953         }
 954         else
 955         {
 956             fprintf(fplog,"Long Range LJ corr.: <C6> %10.4e\n",fr->avcsix[0]);
 957         }
 958     }
 959 }
 960
 961
 962 static void set_bham_b_max(FILE *fplog,t_forcerec *fr,
 963                            const gmx_mtop_t *mtop)
 964 {
 965     const t_atoms *at1,*at2;
 966     int  mt1,mt2,i,j,tpi,tpj,ntypes;
 967     real b,bmin;
 968     real *nbfp;
 969
 970     if (fplog)
 971     {
 972         fprintf(fplog,"Determining largest Buckingham b parameter for table\n");
 973     }
 974     nbfp   = fr->nbfp;
 975     ntypes = fr->ntype;
 976
 977     bmin           = -1;
 978     fr->bham_b_max = 0;
 979     for(mt1=0; mt1<mtop->nmoltype; mt1++)
 980     {
 981         at1 = &mtop->moltype[mt1].atoms;
 982         for(i=0; (i<at1->nr); i++)
 983         {
 984             tpi = at1->atom[i].type;
 985             if (tpi >= ntypes)
 986                 gmx_fatal(FARGS,"Atomtype[%d] = %d, maximum = %d",i,tpi,ntypes);
 987
 988             for(mt2=mt1; mt2<mtop->nmoltype; mt2++)
 989             {
 990                 at2 = &mtop->moltype[mt2].atoms;
 991                 for(j=0; (j<at2->nr); j++) {
 992                     tpj = at2->atom[j].type;
 993                     if (tpj >= ntypes)
 994                     {
 995                         gmx_fatal(FARGS,"Atomtype[%d] = %d, maximum = %d",j,tpj,ntypes);
 996                     }
 997                     b = BHAMB(nbfp,ntypes,tpi,tpj);
 998                     if (b > fr->bham_b_max)
 999                     {
1000                         fr->bham_b_max = b;
1001                     }
1002                     if ((b < bmin) || (bmin==-1))
1003                     {
1004                         bmin = b;
1005                     }
1006                 }
1007             }
1008         }
1009     }
1010     if (fplog)
1011     {
1012         fprintf(fplog,"Buckingham b parameters, min: %g, max: %g\n",
1013                 bmin,fr->bham_b_max);
1014     }
1015 }
1016
1017 static void make_nbf_tables(FILE *fp,const output_env_t oenv,
1018                             t_forcerec *fr,real rtab,
1019                             const t_commrec *cr,
1020                             const char *tabfn,char *eg1,char *eg2,
1021                             t_nblists *nbl)
1022 {
1023   char buf[STRLEN];
1024   int i,j;
1025
1026   if (tabfn == NULL) {
1027     if (debug)
1028       fprintf(debug,"No table file name passed, can not read table, can not do non-bonded interactions\n");
1029     return;
1030   }
1031
1032   sprintf(buf,"%s",tabfn);
1033   if (eg1 && eg2)
1034     /* Append the two energy group names */
1035     sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1,"_%s_%s.%s",
1036             eg1,eg2,ftp2ext(efXVG));
1037   nbl->tab = make_tables(fp,oenv,fr,MASTER(cr),buf,rtab,0);
1038   /* Copy the contents of the table to separate coulomb and LJ tables too,
1039    * to improve cache performance.
1040    */
1041
1042   /* For performance reasons we want
1043    * the table data to be aligned to 16-byte. The pointer could be freed
1044    * but currently isn't.
1045    */
1046   snew_aligned(nbl->vdwtab,8*(nbl->tab.n+1),16);
1047   snew_aligned(nbl->coultab,4*(nbl->tab.n+1),16);
1048
1049   for(i=0; i<=nbl->tab.n; i++) {
1050     for(j=0; j<4; j++)
1051       nbl->coultab[4*i+j] = nbl->tab.tab[12*i+j];
1052     for(j=0; j<8; j++)
1053       nbl->vdwtab [8*i+j] = nbl->tab.tab[12*i+4+j];
1054   }
1055 }
1056
1057 static void count_tables(int ftype1,int ftype2,const gmx_mtop_t *mtop,
1058                          int *ncount,int **count)
1059 {
1060     const gmx_moltype_t *molt;
1061     const t_ilist *il;
1062     int mt,ftype,stride,i,j,tabnr;
1063
1064     for(mt=0; mt<mtop->nmoltype; mt++)
1065     {
1066         molt = &mtop->moltype[mt];
1067         for(ftype=0; ftype<F_NRE; ftype++)
1068         {
1069             if (ftype == ftype1 || ftype == ftype2) {
1070                 il = &molt->ilist[ftype];
1071                 stride = 1 + NRAL(ftype);
1072                 for(i=0; i<il->nr; i+=stride) {
1073                     tabnr = mtop->ffparams.iparams[il->iatoms[i]].tab.table;
1074                     if (tabnr < 0)
1075                         gmx_fatal(FARGS,"A bonded table number is smaller than 0: %d\n",tabnr);
1076                     if (tabnr >= *ncount) {
1077                         srenew(*count,tabnr+1);
1078                         for(j=*ncount; j<tabnr+1; j++)
1079                             (*count)[j] = 0;
1080                         *ncount = tabnr+1;
1081                     }
1082                     (*count)[tabnr]++;
1083                 }
1084             }
1085         }
1086     }
1087 }
1088
1089 static bondedtable_t *make_bonded_tables(FILE *fplog,
1090                                          int ftype1,int ftype2,
1091                                          const gmx_mtop_t *mtop,
1092                                          const char *basefn,const char *tabext)
1093 {
1094     int  i,ncount,*count;
1095     char tabfn[STRLEN];
1096     bondedtable_t *tab;
1097
1098     tab = NULL;
1099
1100     ncount = 0;
1101     count = NULL;
1102     count_tables(ftype1,ftype2,mtop,&ncount,&count);
1103
1104     if (ncount > 0) {
1105         snew(tab,ncount);
1106         for(i=0; i<ncount; i++) {
1107             if (count[i] > 0) {
1108                 sprintf(tabfn,"%s",basefn);
1109                 sprintf(tabfn + strlen(basefn) - strlen(ftp2ext(efXVG)) - 1,"_%s%d.%s",
1110                         tabext,i,ftp2ext(efXVG));
1111                 tab[i] = make_bonded_table(fplog,tabfn,NRAL(ftype1)-2);
1112             }
1113         }
1114         sfree(count);
1115     }
1116
1117     return tab;
1118 }
1119
1120 void forcerec_set_ranges(t_forcerec *fr,
1121                          int ncg_home,int ncg_force,
1122                          int natoms_force,
1123                          int natoms_force_constr,int natoms_f_novirsum)
1124 {
1125     fr->cg0 = 0;
1126     fr->hcg = ncg_home;
1127
1128     /* fr->ncg_force is unused in the standard code,
1129      * but it can be useful for modified code dealing with charge groups.
1130      */
1131     fr->ncg_force           = ncg_force;
1132     fr->natoms_force        = natoms_force;
1133     fr->natoms_force_constr = natoms_force_constr;
1134
1135     if (fr->natoms_force_constr > fr->nalloc_force)
1136     {
1137         fr->nalloc_force = over_alloc_dd(fr->natoms_force_constr);
1138
1139         if (fr->bTwinRange)
1140         {
1141             srenew(fr->f_twin,fr->nalloc_force);
1142         }
1143     }
1144
1145     if (fr->bF_NoVirSum)
1146     {
1147         fr->f_novirsum_n = natoms_f_novirsum;
1148         if (fr->f_novirsum_n > fr->f_novirsum_nalloc)
1149         {
1150             fr->f_novirsum_nalloc = over_alloc_dd(fr->f_novirsum_n);
1151             srenew(fr->f_novirsum_alloc,fr->f_novirsum_nalloc);
1152         }
1153     }
1154     else
1155     {
1156         fr->f_novirsum_n = 0;
1157     }
1158 }
1159
1160 static real cutoff_inf(real cutoff)
1161 {
1162     if (cutoff == 0)
1163     {
1164         cutoff = GMX_CUTOFF_INF;
1165     }
1166
1167     return cutoff;
1168 }
1169
1170 static void make_adress_tf_tables(FILE *fp,const output_env_t oenv,
1171                             t_forcerec *fr,const t_inputrec *ir,
1172                             const char *tabfn, const gmx_mtop_t *mtop,
1173                             matrix     box)
1174 {
1175   char buf[STRLEN];
1176   int i,j;
1177
1178   if (tabfn == NULL) {
1179         gmx_fatal(FARGS,"No thermoforce table file given. Use -tabletf to specify a file\n");
1180     return;
1181   }
1182
1183   snew(fr->atf_tabs, ir->adress->n_tf_grps);
1184
1185   for (i=0; i<ir->adress->n_tf_grps; i++){
1186     j = ir->adress->tf_table_index[i]; /* get energy group index */
1187     sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1,"tf_%s.%s",
1188         *(mtop->groups.grpname[mtop->groups.grps[egcENER].nm_ind[j]]) ,ftp2ext(efXVG));
1189     printf("loading tf table for energygrp index %d from %s\n", ir->adress->tf_table_index[j], buf);
1190     fr->atf_tabs[i] = make_atf_table(fp,oenv,fr,buf, box);
1191   }
1192
1193 }
1194
1195 gmx_bool can_use_allvsall(const t_inputrec *ir, const gmx_mtop_t *mtop,
1196                       gmx_bool bPrintNote,t_commrec *cr,FILE *fp)
1197 {
1198     gmx_bool bAllvsAll;
1199
1200     bAllvsAll =
1201         (
1202          ir->rlist==0            &&
1203          ir->rcoulomb==0         &&
1204          ir->rvdw==0             &&
1205          ir->ePBC==epbcNONE      &&
1206          ir->vdwtype==evdwCUT    &&
1207          ir->coulombtype==eelCUT &&
1208          ir->efep==efepNO        &&
1209          (ir->implicit_solvent == eisNO ||
1210           (ir->implicit_solvent==eisGBSA && (ir->gb_algorithm==egbSTILL ||
1211                                              ir->gb_algorithm==egbHCT   ||
1212                                              ir->gb_algorithm==egbOBC))) &&
1213          getenv("GMX_NO_ALLVSALL") == NULL
1214             );
1215
1216     if (bAllvsAll && ir->opts.ngener > 1)
1217     {
1218         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";
1219
1220         if (bPrintNote)
1221         {
1222             if (MASTER(cr))
1223             {
1224                 fprintf(stderr,"\n%s\n",note);
1225             }
1226             if (fp != NULL)
1227             {
1228                 fprintf(fp,"\n%s\n",note);
1229             }
1230         }
1231         bAllvsAll = FALSE;
1232     }
1233
1234     if(bAllvsAll && fp && MASTER(cr))
1235     {
1236         fprintf(fp,"\nUsing accelerated all-vs-all kernels.\n\n");
1237     }
1238
1239     return bAllvsAll;
1240 }
1241
1242
1243
1244 void init_forcerec(FILE *fp,
1245                    const output_env_t oenv,
1246                    t_forcerec *fr,
1247                    t_fcdata   *fcd,
1248                    const t_inputrec *ir,
1249                    const gmx_mtop_t *mtop,
1250                    const t_commrec  *cr,
1251                    matrix     box,
1252                    gmx_bool       bMolEpot,
1253                    const char *tabfn,
1254                    const char *tabafn,
1255                    const char *tabpfn,
1256                    const char *tabbfn,
1257                    gmx_bool       bNoSolvOpt,
1258                    real       print_force)
1259 {
1260     int     i,j,m,natoms,ngrp,negp_pp,negptable,egi,egj;
1261     real    rtab;
1262     char    *env;
1263     double  dbl;
1264     rvec    box_size;
1265     const t_block *cgs;
1266     gmx_bool    bGenericKernelOnly;
1267     gmx_bool    bTab,bSep14tab,bNormalnblists;
1268     t_nblists *nbl;
1269     int     *nm_ind,egp_flags;
1270
1271     gmx_detectcpu(&fr->cpu_information);
1272     if(MASTER(cr))
1273     {
1274         /* Only print warnings from master */
1275         gmx_detectcpu_check_acceleration(fr->cpu_information,fp);
1276     }
1277
1278     /* By default we turn acceleration on, but it might be turned off further down... */
1279     fr->use_acceleration = TRUE;
1280
1281     fr->bDomDec = DOMAINDECOMP(cr);
1282
1283     natoms = mtop->natoms;
1284
1285     if (check_box(ir->ePBC,box))
1286     {
1287         gmx_fatal(FARGS,check_box(ir->ePBC,box));
1288     }
1289
1290     /* Test particle insertion ? */
1291     if (EI_TPI(ir->eI)) {
1292         /* Set to the size of the molecule to be inserted (the last one) */
1293         /* Because of old style topologies, we have to use the last cg
1294          * instead of the last molecule type.
1295          */
1296         cgs = &mtop->moltype[mtop->molblock[mtop->nmolblock-1].type].cgs;
1297         fr->n_tpi = cgs->index[cgs->nr] - cgs->index[cgs->nr-1];
1298         if (fr->n_tpi != mtop->mols.index[mtop->mols.nr] - mtop->mols.index[mtop->mols.nr-1]) {
1299             gmx_fatal(FARGS,"The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group.");
1300         }
1301     } else {
1302         fr->n_tpi = 0;
1303     }
1304
1305     /* Copy AdResS parameters */
1306     if (ir->bAdress) {
1307       fr->adress_type     = ir->adress->type;
1308       fr->adress_const_wf = ir->adress->const_wf;
1309       fr->adress_ex_width = ir->adress->ex_width;
1310       fr->adress_hy_width = ir->adress->hy_width;
1311       fr->adress_icor     = ir->adress->icor;
1312       fr->adress_site     = ir->adress->site;
1313       fr->adress_ex_forcecap = ir->adress->ex_forcecap;
1314       fr->adress_do_hybridpairs = ir->adress->do_hybridpairs;
1315
1316
1317       snew(fr->adress_group_explicit , ir->adress->n_energy_grps);
1318       for (i=0; i< ir->adress->n_energy_grps; i++){
1319           fr->adress_group_explicit[i]= ir->adress->group_explicit[i];
1320       }
1321
1322       fr->n_adress_tf_grps = ir->adress->n_tf_grps;
1323       snew(fr->adress_tf_table_index, fr->n_adress_tf_grps);
1324       for (i=0; i< fr->n_adress_tf_grps; i++){
1325           fr->adress_tf_table_index[i]= ir->adress->tf_table_index[i];
1326       }
1327       copy_rvec(ir->adress->refs,fr->adress_refs);
1328     } else {
1329       fr->adress_type = eAdressOff;
1330       fr->adress_do_hybridpairs = FALSE;
1331     }
1332
1333     /* Copy the user determined parameters */
1334     fr->userint1 = ir->userint1;
1335     fr->userint2 = ir->userint2;
1336     fr->userint3 = ir->userint3;
1337     fr->userint4 = ir->userint4;
1338     fr->userreal1 = ir->userreal1;
1339     fr->userreal2 = ir->userreal2;
1340     fr->userreal3 = ir->userreal3;
1341     fr->userreal4 = ir->userreal4;
1342
1343     /* Shell stuff */
1344     fr->fc_stepsize = ir->fc_stepsize;
1345
1346     /* Free energy */
1347     fr->efep       = ir->efep;
1348     fr->sc_alphavdw = ir->fepvals->sc_alpha;
1349     if (ir->fepvals->bScCoul)
1350     {
1351         fr->sc_alphacoul = ir->fepvals->sc_alpha;
1352         fr->sc_sigma6_min = pow(ir->fepvals->sc_sigma_min,6);
1353     }
1354     else
1355     {
1356         fr->sc_alphacoul = 0;
1357         fr->sc_sigma6_min = 0; /* only needed when bScCoul is on */
1358     }
1359     fr->sc_power   = ir->fepvals->sc_power;
1360     fr->sc_r_power   = ir->fepvals->sc_r_power;
1361     fr->sc_sigma6_def = pow(ir->fepvals->sc_sigma,6);
1362
1363     env = getenv("GMX_SCSIGMA_MIN");
1364     if (env != NULL)
1365     {
1366         dbl = 0;
1367         sscanf(env,"%lf",&dbl);
1368         fr->sc_sigma6_min = pow(dbl,6);
1369         if (fp)
1370         {
1371             fprintf(fp,"Setting the minimum soft core sigma to %g nm\n",dbl);
1372         }
1373     }
1374
1375     bGenericKernelOnly = FALSE;
1376     if (getenv("GMX_NB_GENERIC") != NULL)
1377     {
1378         if (fp != NULL)
1379         {
1380             fprintf(fp,
1381                     "Found environment variable GMX_NB_GENERIC.\n"
1382                     "Disabling interaction-specific nonbonded kernels.\n\n");
1383         }
1384         bGenericKernelOnly = TRUE;
1385         bNoSolvOpt         = TRUE;
1386     }
1387
1388     if( (getenv("GMX_DISABLE_ACCELERATION") != NULL) || (getenv("GMX_NOOPTIMIZEDKERNELS") != NULL) )
1389     {
1390         fr->use_acceleration = FALSE;
1391         if (fp != NULL)
1392         {
1393             fprintf(fp,
1394                     "\nFound environment variable GMX_DISABLE_ACCELERATION.\n"
1395                     "Disabling all architecture-specific (e.g. SSE2/SSE4/AVX) routines.\n\n");
1396         }
1397     }
1398
1399     /* Check if we can/should do all-vs-all kernels */
1400     fr->bAllvsAll       = can_use_allvsall(ir,mtop,FALSE,NULL,NULL);
1401     fr->AllvsAll_work   = NULL;
1402     fr->AllvsAll_workgb = NULL;
1403
1404
1405
1406     /* Neighbour searching stuff */
1407     fr->bGrid      = (ir->ns_type == ensGRID);
1408     fr->ePBC       = ir->ePBC;
1409     fr->bMolPBC    = ir->bPeriodicMols;
1410     fr->rc_scaling = ir->refcoord_scaling;
1411     copy_rvec(ir->posres_com,fr->posres_com);
1412     copy_rvec(ir->posres_comB,fr->posres_comB);
1413     fr->rlist      = cutoff_inf(ir->rlist);
1414     fr->rlistlong  = cutoff_inf(ir->rlistlong);
1415     fr->eeltype    = ir->coulombtype;
1416     fr->vdwtype    = ir->vdwtype;
1417
1418     fr->bTwinRange = fr->rlistlong > fr->rlist;
1419     fr->bEwald     = (EEL_PME(fr->eeltype) || fr->eeltype==eelEWALD);
1420
1421     fr->reppow     = mtop->ffparams.reppow;
1422     fr->bvdwtab    = (fr->vdwtype != evdwCUT ||
1423                       !gmx_within_tol(fr->reppow,12.0,10*GMX_DOUBLE_EPS));
1424     fr->bcoultab   = (!(fr->eeltype == eelCUT || EEL_RF(fr->eeltype)) ||
1425                       fr->eeltype == eelRF_ZERO);
1426
1427     if (getenv("GMX_FORCE_TABLES"))
1428     {
1429         fr->bvdwtab  = TRUE;
1430         fr->bcoultab = TRUE;
1431     }
1432
1433     if (fp) {
1434         fprintf(fp,"Table routines are used for coulomb: %s\n",bool_names[fr->bcoultab]);
1435         fprintf(fp,"Table routines are used for vdw:     %s\n",bool_names[fr->bvdwtab ]);
1436     }
1437
1438     /* Tables are used for direct ewald sum */
1439     if(fr->bEwald)
1440     {
1441         if (EEL_PME(ir->coulombtype))
1442         {
1443             if (fp)
1444                 fprintf(fp,"Will do PME sum in reciprocal space.\n");
1445             if (ir->coulombtype == eelP3M_AD)
1446             {
1447                 please_cite(fp,"Hockney1988");
1448                 please_cite(fp,"Ballenegger2012");
1449             }
1450             else
1451             {
1452                 please_cite(fp,"Essmann95a");
1453             }
1454
1455             if (ir->ewald_geometry == eewg3DC)
1456             {
1457                 if (fp)
1458                 {
1459                     fprintf(fp,"Using the Ewald3DC correction for systems with a slab geometry.\n");
1460                 }
1461                 please_cite(fp,"In-Chul99a");
1462             }
1463         }
1464         fr->ewaldcoeff=calc_ewaldcoeff(ir->rcoulomb, ir->ewald_rtol);
1465         init_ewald_tab(&(fr->ewald_table), cr, ir, fp);
1466         if (fp)
1467         {
1468             fprintf(fp,"Using a Gaussian width (1/beta) of %g nm for Ewald\n",
1469                     1/fr->ewaldcoeff);
1470         }
1471     }
1472
1473     /* Electrostatics */
1474     fr->epsilon_r  = ir->epsilon_r;
1475     fr->epsilon_rf = ir->epsilon_rf;
1476     fr->fudgeQQ    = mtop->ffparams.fudgeQQ;
1477     fr->rcoulomb_switch = ir->rcoulomb_switch;
1478     fr->rcoulomb        = cutoff_inf(ir->rcoulomb);
1479
1480     /* Parameters for generalized RF */
1481     fr->zsquare = 0.0;
1482     fr->temp    = 0.0;
1483
1484     if (fr->eeltype == eelGRF)
1485     {
1486         init_generalized_rf(fp,mtop,ir,fr);
1487     }
1488     else if (fr->eeltype == eelSHIFT)
1489     {
1490         for(m=0; (m<DIM); m++)
1491             box_size[m]=box[m][m];
1492
1493         if ((fr->eeltype == eelSHIFT && fr->rcoulomb > fr->rcoulomb_switch))
1494             set_shift_consts(fp,fr->rcoulomb_switch,fr->rcoulomb,box_size,fr);
1495     }
1496
1497     fr->bF_NoVirSum = (EEL_FULL(fr->eeltype) ||
1498                        gmx_mtop_ftype_count(mtop,F_POSRES) > 0 ||
1499                        IR_ELEC_FIELD(*ir) ||
1500                        (fr->adress_icor != eAdressICOff)
1501                       );
1502
1503     /* Mask that says whether or not this NBF list should be computed */
1504     /*  if (fr->bMask == NULL) {
1505         ngrp = ir->opts.ngener*ir->opts.ngener;
1506         snew(fr->bMask,ngrp);*/
1507     /* Defaults to always */
1508     /*    for(i=0; (i<ngrp); i++)
1509           fr->bMask[i] = TRUE;
1510           }*/
1511
1512     if (ncg_mtop(mtop) > fr->cg_nalloc && !DOMAINDECOMP(cr)) {
1513         /* Count the total number of charge groups */
1514         fr->cg_nalloc = ncg_mtop(mtop);
1515         srenew(fr->cg_cm,fr->cg_nalloc);
1516     }
1517     if (fr->shift_vec == NULL)
1518         snew(fr->shift_vec,SHIFTS);
1519
1520     if (fr->fshift == NULL)
1521         snew(fr->fshift,SHIFTS);
1522
1523     if (fr->nbfp == NULL) {
1524         fr->ntype = mtop->ffparams.atnr;
1525         fr->bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
1526         fr->nbfp  = mk_nbfp(&mtop->ffparams,fr->bBHAM);
1527     }
1528
1529     /* Copy the energy group exclusions */
1530     fr->egp_flags = ir->opts.egp_flags;
1531
1532     /* Van der Waals stuff */
1533     fr->rvdw        = cutoff_inf(ir->rvdw);
1534     fr->rvdw_switch = ir->rvdw_switch;
1535     if ((fr->vdwtype != evdwCUT) && (fr->vdwtype != evdwUSER) && !fr->bBHAM) {
1536         if (fr->rvdw_switch >= fr->rvdw)
1537             gmx_fatal(FARGS,"rvdw_switch (%f) must be < rvdw (%f)",
1538                       fr->rvdw_switch,fr->rvdw);
1539         if (fp)
1540             fprintf(fp,"Using %s Lennard-Jones, switch between %g and %g nm\n",
1541                     (fr->eeltype==eelSWITCH) ? "switched":"shifted",
1542                     fr->rvdw_switch,fr->rvdw);
1543     }
1544
1545     if (fr->bBHAM && (fr->vdwtype == evdwSHIFT || fr->vdwtype == evdwSWITCH))
1546         gmx_fatal(FARGS,"Switch/shift interaction not supported with Buckingham");
1547
1548     if (fp)
1549         fprintf(fp,"Cut-off's:   NS: %g   Coulomb: %g   %s: %g\n",
1550                 fr->rlist,fr->rcoulomb,fr->bBHAM ? "BHAM":"LJ",fr->rvdw);
1551
1552     fr->eDispCorr = ir->eDispCorr;
1553     if (ir->eDispCorr != edispcNO)
1554     {
1555         set_avcsixtwelve(fp,fr,mtop);
1556     }
1557
1558     if (fr->bBHAM)
1559     {
1560         set_bham_b_max(fp,fr,mtop);
1561     }
1562
1563     fr->bGB = (ir->implicit_solvent == eisGBSA);
1564         fr->gb_epsilon_solvent = ir->gb_epsilon_solvent;
1565
1566     /* Copy the GBSA data (radius, volume and surftens for each
1567      * atomtype) from the topology atomtype section to forcerec.
1568      */
1569     snew(fr->atype_radius,fr->ntype);
1570     snew(fr->atype_vol,fr->ntype);
1571     snew(fr->atype_surftens,fr->ntype);
1572     snew(fr->atype_gb_radius,fr->ntype);
1573     snew(fr->atype_S_hct,fr->ntype);
1574
1575     if (mtop->atomtypes.nr > 0)
1576     {
1577         for(i=0;i<fr->ntype;i++)
1578             fr->atype_radius[i] =mtop->atomtypes.radius[i];
1579         for(i=0;i<fr->ntype;i++)
1580             fr->atype_vol[i] = mtop->atomtypes.vol[i];
1581         for(i=0;i<fr->ntype;i++)
1582             fr->atype_surftens[i] = mtop->atomtypes.surftens[i];
1583         for(i=0;i<fr->ntype;i++)
1584             fr->atype_gb_radius[i] = mtop->atomtypes.gb_radius[i];
1585         for(i=0;i<fr->ntype;i++)
1586             fr->atype_S_hct[i] = mtop->atomtypes.S_hct[i];
1587     }
1588
1589         /* Generate the GB table if needed */
1590         if(fr->bGB)
1591         {
1592 #ifdef GMX_DOUBLE
1593                 fr->gbtabscale=2000;
1594 #else
1595                 fr->gbtabscale=500;
1596 #endif
1597
1598                 fr->gbtabr=100;
1599                 fr->gbtab=make_gb_table(fp,oenv,fr,tabpfn,fr->gbtabscale);
1600
1601         init_gb(&fr->born,cr,fr,ir,mtop,ir->rgbradii,ir->gb_algorithm);
1602
1603         /* Copy local gb data (for dd, this is done in dd_partition_system) */
1604         if (!DOMAINDECOMP(cr))
1605         {
1606             make_local_gb(cr,fr->born,ir->gb_algorithm);
1607         }
1608     }
1609
1610     /* Set the charge scaling */
1611     if (fr->epsilon_r != 0)
1612         fr->epsfac = ONE_4PI_EPS0/fr->epsilon_r;
1613     else
1614         /* eps = 0 is infinite dieletric: no coulomb interactions */
1615         fr->epsfac = 0;
1616
1617     /* Reaction field constants */
1618     if (EEL_RF(fr->eeltype))
1619         calc_rffac(fp,fr->eeltype,fr->epsilon_r,fr->epsilon_rf,
1620                    fr->rcoulomb,fr->temp,fr->zsquare,box,
1621                    &fr->kappa,&fr->k_rf,&fr->c_rf);
1622
1623     set_chargesum(fp,fr,mtop);
1624
1625     /* if we are using LR electrostatics, and they are tabulated,
1626      * the tables will contain modified coulomb interactions.
1627      * Since we want to use the non-shifted ones for 1-4
1628      * coulombic interactions, we must have an extra set of tables.
1629      */
1630
1631     /* Construct tables.
1632      * A little unnecessary to make both vdw and coul tables sometimes,
1633      * but what the heck... */
1634
1635     bTab = fr->bcoultab || fr->bvdwtab;
1636
1637     bSep14tab = ((!bTab || fr->eeltype!=eelCUT || fr->vdwtype!=evdwCUT ||
1638                   fr->bBHAM) &&
1639                  (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
1640                   gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0 ||
1641                   gmx_mtop_ftype_count(mtop,F_LJC_PAIRS_NB) > 0));
1642
1643     negp_pp = ir->opts.ngener - ir->nwall;
1644     negptable = 0;
1645     if (!bTab) {
1646         bNormalnblists = TRUE;
1647         fr->nnblists = 1;
1648     } else {
1649         bNormalnblists = (ir->eDispCorr != edispcNO);
1650         for(egi=0; egi<negp_pp; egi++) {
1651             for(egj=egi;  egj<negp_pp; egj++) {
1652                 egp_flags = ir->opts.egp_flags[GID(egi,egj,ir->opts.ngener)];
1653                 if (!(egp_flags & EGP_EXCL)) {
1654                     if (egp_flags & EGP_TABLE) {
1655                         negptable++;
1656                     } else {
1657                         bNormalnblists = TRUE;
1658                     }
1659                 }
1660             }
1661         }
1662         if (bNormalnblists) {
1663             fr->nnblists = negptable + 1;
1664         } else {
1665             fr->nnblists = negptable;
1666         }
1667         if (fr->nnblists > 1)
1668             snew(fr->gid2nblists,ir->opts.ngener*ir->opts.ngener);
1669     }
1670     snew(fr->nblists,fr->nnblists);
1671
1672     /* This code automatically gives table length tabext without cut-off's,
1673      * in that case grompp should already have checked that we do not need
1674      * normal tables and we only generate tables for 1-4 interactions.
1675      */
1676     rtab = ir->rlistlong + ir->tabext;
1677
1678     if (bTab) {
1679         /* make tables for ordinary interactions */
1680         if (bNormalnblists) {
1681             make_nbf_tables(fp,oenv,fr,rtab,cr,tabfn,NULL,NULL,&fr->nblists[0]);
1682             if (!bSep14tab)
1683                 fr->tab14 = fr->nblists[0].tab;
1684             m = 1;
1685         } else {
1686             m = 0;
1687         }
1688         if (negptable > 0) {
1689             /* Read the special tables for certain energy group pairs */
1690             nm_ind = mtop->groups.grps[egcENER].nm_ind;
1691             for(egi=0; egi<negp_pp; egi++) {
1692                 for(egj=egi;  egj<negp_pp; egj++) {
1693                     egp_flags = ir->opts.egp_flags[GID(egi,egj,ir->opts.ngener)];
1694                     if ((egp_flags & EGP_TABLE) && !(egp_flags & EGP_EXCL)) {
1695                         nbl = &(fr->nblists[m]);
1696                         if (fr->nnblists > 1) {
1697                             fr->gid2nblists[GID(egi,egj,ir->opts.ngener)] = m;
1698                         }
1699                         /* Read the table file with the two energy groups names appended */
1700                         make_nbf_tables(fp,oenv,fr,rtab,cr,tabfn,
1701                                         *mtop->groups.grpname[nm_ind[egi]],
1702                                         *mtop->groups.grpname[nm_ind[egj]],
1703                                         &fr->nblists[m]);
1704                         m++;
1705                     } else if (fr->nnblists > 1) {
1706                         fr->gid2nblists[GID(egi,egj,ir->opts.ngener)] = 0;
1707                     }
1708                 }
1709             }
1710         }
1711     }
1712     if (bSep14tab)
1713     {
1714         /* generate extra tables with plain Coulomb for 1-4 interactions only */
1715         fr->tab14 = make_tables(fp,oenv,fr,MASTER(cr),tabpfn,rtab,
1716                                 GMX_MAKETABLES_14ONLY);
1717     }
1718
1719     /* Read AdResS Thermo Force table if needed */
1720     if(fr->adress_icor == eAdressICThermoForce)
1721     {
1722         /* old todo replace */
1723
1724         if (ir->adress->n_tf_grps > 0){
1725             make_adress_tf_tables(fp,oenv,fr,ir,tabfn, mtop, box);
1726
1727         }else{
1728             /* load the default table */
1729             snew(fr->atf_tabs, 1);
1730             fr->atf_tabs[DEFAULT_TF_TABLE] = make_atf_table(fp,oenv,fr,tabafn, box);
1731         }
1732     }
1733
1734     /* Wall stuff */
1735     fr->nwall = ir->nwall;
1736     if (ir->nwall && ir->wall_type==ewtTABLE)
1737     {
1738         make_wall_tables(fp,oenv,ir,tabfn,&mtop->groups,fr);
1739     }
1740
1741     if (fcd && tabbfn) {
1742         fcd->bondtab  = make_bonded_tables(fp,
1743                                            F_TABBONDS,F_TABBONDSNC,
1744                                            mtop,tabbfn,"b");
1745         fcd->angletab = make_bonded_tables(fp,
1746                                            F_TABANGLES,-1,
1747                                            mtop,tabbfn,"a");
1748         fcd->dihtab   = make_bonded_tables(fp,
1749                                            F_TABDIHS,-1,
1750                                            mtop,tabbfn,"d");
1751     } else {
1752         if (debug)
1753             fprintf(debug,"No fcdata or table file name passed, can not read table, can not do bonded interactions\n");
1754     }
1755
1756     if (ir->eDispCorr != edispcNO)
1757     {
1758         calc_enervirdiff(fp,ir->eDispCorr,fr);
1759     }
1760
1761     /* QM/MM initialization if requested
1762      */
1763     if (ir->bQMMM)
1764     {
1765         fprintf(stderr,"QM/MM calculation requested.\n");
1766     }
1767
1768     fr->bQMMM      = ir->bQMMM;
1769     fr->qr         = mk_QMMMrec();
1770
1771     /* Set all the static charge group info */
1772     fr->cginfo_mb = init_cginfo_mb(fp,mtop,fr,bNoSolvOpt);
1773
1774     if (DOMAINDECOMP(cr)) {
1775         fr->cginfo = NULL;
1776     } else {
1777         fr->cginfo = cginfo_expand(mtop->nmolblock,fr->cginfo_mb);
1778     }
1779
1780     if (!DOMAINDECOMP(cr))
1781     {
1782         /* When using particle decomposition, the effect of the second argument,
1783          * which sets fr->hcg, is corrected later in do_md and init_em.
1784          */
1785         forcerec_set_ranges(fr,ncg_mtop(mtop),ncg_mtop(mtop),
1786                             mtop->natoms,mtop->natoms,mtop->natoms);
1787     }
1788
1789     fr->print_force = print_force;
1790
1791
1792     /* coarse load balancing vars */
1793     fr->t_fnbf=0.;
1794     fr->t_wait=0.;
1795     fr->timesteps=0;
1796
1797     /* Initialize neighbor search */
1798     init_ns(fp,cr,&fr->ns,fr,mtop,box);
1799
1800     if (cr->duty & DUTY_PP){
1801         gmx_setup_kernels(fp,fr,bGenericKernelOnly);
1802         if (ir->bAdress)
1803             gmx_setup_adress_kernels(fp,bGenericKernelOnly);
1804     }
1805 }
1806
1807 #define pr_real(fp,r) fprintf(fp,"%s: %e\n",#r,r)
1808 #define pr_int(fp,i)  fprintf((fp),"%s: %d\n",#i,i)
1809 #define pr_bool(fp,b) fprintf((fp),"%s: %s\n",#b,bool_names[b])
1810
1811 void pr_forcerec(FILE *fp,t_forcerec *fr,t_commrec *cr)
1812 {
1813   int i;
1814
1815   pr_real(fp,fr->rlist);
1816   pr_real(fp,fr->rcoulomb);
1817   pr_real(fp,fr->fudgeQQ);
1818   pr_bool(fp,fr->bGrid);
1819   pr_bool(fp,fr->bTwinRange);
1820   /*pr_int(fp,fr->cg0);
1821     pr_int(fp,fr->hcg);*/
1822   for(i=0; i<fr->nnblists; i++)
1823     pr_int(fp,fr->nblists[i].tab.n);
1824   pr_real(fp,fr->rcoulomb_switch);
1825   pr_real(fp,fr->rcoulomb);
1826
1827   fflush(fp);
1828 }