src/gmxlib/nonbonded/nonbonded.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  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #ifdef GMX_THREAD_SHM_FDECOMP
  41 #include <thread_mpi.h>
  42 #endif
  43
  44
  45 #include <stdio.h>
  46 #include <stdlib.h>
  47 #include "typedefs.h"
  48 #include "txtdump.h"
  49 #include "smalloc.h"
  50 #include "ns.h"
  51 #include "vec.h"
  52 #include "maths.h"
  53 #include "macros.h"
  54 #include "force.h"
  55 #include "names.h"
  56 #include "main.h"
  57 #include "xvgr.h"
  58 #include "gmx_fatal.h"
  59 #include "physics.h"
  60 #include "force.h"
  61 #include "bondf.h"
  62 #include "nrnb.h"
  63 #include "smalloc.h"
  64 #include "nonbonded.h"
  65
  66 #include "nb_kernel_c/nb_kernel_c.h"
  67 #include "nb_kernel_adress_c/nb_kernel_c_adress.h"
  68 #include "nb_free_energy.h"
  69 #include "nb_generic.h"
  70 #include "nb_generic_cg.h"
  71 #include "nb_generic_adress.h"
  72
  73
  74 /* 1,4 interactions uses kernel 330 directly */
  75 #include "nb_kernel_c/nb_kernel330.h"
  76 #include "nb_kernel_adress_c/nb_kernel330_adress.h"
  77
  78 #ifdef GMX_PPC_ALTIVEC
  79 #include "nb_kernel_ppc_altivec/nb_kernel_ppc_altivec.h"
  80 #endif
  81
  82 #if defined(GMX_IA32_SSE)
  83 #include "nb_kernel_ia32_sse/nb_kernel_ia32_sse.h"
  84 #endif
  85
  86 #if defined(GMX_IA32_SSE2)
  87 #include "nb_kernel_ia32_sse2/nb_kernel_ia32_sse2.h"
  88 #endif
  89
  90 #if defined(GMX_X86_64_SSE)
  91 #include "nb_kernel_x86_64_sse/nb_kernel_x86_64_sse.h"
  92 #endif
  93
  94 #if defined(GMX_X86_64_SSE2)
  95 #include "nb_kernel_x86_64_sse2/nb_kernel_x86_64_sse2.h"
  96 #endif
  97
  98 #if defined(GMX_SSE2)
  99 #  ifdef GMX_DOUBLE
 100 #    include "nb_kernel_sse2_double/nb_kernel_sse2_double.h"
 101 #  else
 102 #    include "nb_kernel_sse2_single/nb_kernel_sse2_single.h"
 103 #  endif
 104 #endif
 105
 106 #if defined(GMX_FORTRAN)
 107 #  ifdef GMX_DOUBLE
 108 #    include "nb_kernel_f77_double/nb_kernel_f77_double.h"
 109 #  else
 110 #    include "nb_kernel_f77_single/nb_kernel_f77_single.h"
 111 #  endif
 112 #endif
 113
 114 #if (defined GMX_IA64_ASM && defined GMX_DOUBLE)
 115 #include "nb_kernel_ia64_double/nb_kernel_ia64_double.h"
 116 #endif
 117
 118 #if (defined GMX_IA64_ASM && !defined GMX_DOUBLE)
 119 #include "nb_kernel_ia64_single/nb_kernel_ia64_single.h"
 120 #endif
 121
 122 #ifdef GMX_POWER6
 123 #include "nb_kernel_power6/nb_kernel_power6.h"
 124 #endif
 125
 126 #ifdef GMX_BLUEGENE
 127 #include "nb_kernel_bluegene/nb_kernel_bluegene.h"
 128 #endif
 129
 130
 131
 132 enum { TABLE_NONE, TABLE_COMBINED, TABLE_COUL, TABLE_VDW, TABLE_NR };
 133
 134
 135 /* Table version for each kernel.
 136  */
 137 static const int
 138 nb_kernel_table[eNR_NBKERNEL_NR] =
 139 {
 140   TABLE_NONE,     /* kernel010 */
 141   TABLE_NONE,     /* kernel020 */
 142   TABLE_VDW,      /* kernel030 */
 143   TABLE_NONE,     /* kernel100 */
 144   TABLE_NONE,     /* kernel101 */
 145   TABLE_NONE,     /* kernel102 */
 146   TABLE_NONE,     /* kernel103 */
 147   TABLE_NONE,     /* kernel104 */
 148   TABLE_NONE,     /* kernel110 */
 149   TABLE_NONE,     /* kernel111 */
 150   TABLE_NONE,     /* kernel112 */
 151   TABLE_NONE,     /* kernel113 */
 152   TABLE_NONE,     /* kernel114 */
 153   TABLE_NONE,     /* kernel120 */
 154   TABLE_NONE,     /* kernel121 */
 155   TABLE_NONE,     /* kernel122 */
 156   TABLE_NONE,     /* kernel123 */
 157   TABLE_NONE,     /* kernel124 */
 158   TABLE_VDW,      /* kernel130 */
 159   TABLE_VDW,      /* kernel131 */
 160   TABLE_VDW,      /* kernel132 */
 161   TABLE_VDW,      /* kernel133 */
 162   TABLE_VDW,      /* kernel134 */
 163   TABLE_NONE,     /* kernel200 */
 164   TABLE_NONE,     /* kernel201 */
 165   TABLE_NONE,     /* kernel202 */
 166   TABLE_NONE,     /* kernel203 */
 167   TABLE_NONE,     /* kernel204 */
 168   TABLE_NONE,     /* kernel210 */
 169   TABLE_NONE,     /* kernel211 */
 170   TABLE_NONE,     /* kernel212 */
 171   TABLE_NONE,     /* kernel213 */
 172   TABLE_NONE,     /* kernel214 */
 173   TABLE_NONE,     /* kernel220 */
 174   TABLE_NONE,     /* kernel221 */
 175   TABLE_NONE,     /* kernel222 */
 176   TABLE_NONE,     /* kernel223 */
 177   TABLE_NONE,     /* kernel224 */
 178   TABLE_VDW,      /* kernel230 */
 179   TABLE_VDW,      /* kernel231 */
 180   TABLE_VDW,      /* kernel232 */
 181   TABLE_VDW,      /* kernel233 */
 182   TABLE_VDW,      /* kernel234 */
 183   TABLE_COUL,     /* kernel300 */
 184   TABLE_COUL,     /* kernel301 */
 185   TABLE_COUL,     /* kernel302 */
 186   TABLE_COUL,     /* kernel303 */
 187   TABLE_COUL,     /* kernel304 */
 188   TABLE_COUL,     /* kernel310 */
 189   TABLE_COUL,     /* kernel311 */
 190   TABLE_COUL,     /* kernel312 */
 191   TABLE_COUL,     /* kernel313 */
 192   TABLE_COUL,     /* kernel314 */
 193   TABLE_COUL,     /* kernel320 */
 194   TABLE_COUL,     /* kernel321 */
 195   TABLE_COUL,     /* kernel322 */
 196   TABLE_COUL,     /* kernel323 */
 197   TABLE_COUL,     /* kernel324 */
 198   TABLE_COMBINED, /* kernel330 */
 199   TABLE_COMBINED, /* kernel331 */
 200   TABLE_COMBINED, /* kernel332 */
 201   TABLE_COMBINED, /* kernel333 */
 202   TABLE_COMBINED, /* kernel334 */
 203   TABLE_NONE,     /* kernel400 */
 204   TABLE_NONE,     /* kernel410 */
 205   TABLE_VDW       /* kernel430 */
 206 };
 207
 208
 209
 210 static nb_kernel_t **
 211 nb_kernel_list = NULL;
 212
 213 static nb_adress_kernel_t **
 214 nb_kernel_list_adress = NULL;
 215
 216 void
 217 gmx_setup_kernels(FILE *fplog,gmx_bool bGenericKernelOnly)
 218 {
 219     int i;
 220
 221     snew(nb_kernel_list,eNR_NBKERNEL_NR);
 222
 223     /* Note that later calls overwrite earlier, so the preferred (fastest)
 224      * version should be at the end. For instance, we call SSE after 3DNow.
 225      */
 226
 227     for(i=0; i<eNR_NBKERNEL_NR; i++)
 228     {
 229         nb_kernel_list[i] = NULL;
 230     }
 231
 232     if (bGenericKernelOnly)
 233     {
 234         return;
 235     }
 236
 237         if(fplog)
 238     {
 239             fprintf(fplog,"Configuring nonbonded kernels...\n");
 240     }
 241
 242     nb_kernel_setup(fplog,nb_kernel_list);
 243
 244     if(getenv("GMX_NOOPTIMIZEDKERNELS") != NULL)
 245     {
 246         return;
 247     }
 248
 249     /* Setup kernels. The last called setup routine will overwrite earlier assignments,
 250          * so we should e.g. test SSE3 support _after_ SSE2 support,
 251      * and call e.g. Fortran setup before SSE.
 252          */
 253
 254 #if defined(GMX_FORTRAN) && defined(GMX_DOUBLE)
 255     nb_kernel_setup_f77_double(fplog,nb_kernel_list);
 256 #endif
 257
 258 #if defined(GMX_FORTRAN) && !defined(GMX_DOUBLE)
 259     nb_kernel_setup_f77_single(fplog,nb_kernel_list);
 260 #endif
 261
 262 #ifdef GMX_BLUEGENE
 263     nb_kernel_setup_bluegene(fplog,nb_kernel_list);
 264 #endif
 265
 266 #ifdef GMX_POWER6
 267     nb_kernel_setup_power6(fplog,nb_kernel_list);
 268 #endif
 269
 270 #ifdef GMX_PPC_ALTIVEC
 271     nb_kernel_setup_ppc_altivec(fplog,nb_kernel_list);
 272 #endif
 273
 274 #if defined(GMX_IA32_SSE)
 275     nb_kernel_setup_ia32_sse(fplog,nb_kernel_list);
 276 #endif
 277
 278 #if defined(GMX_IA32_SSE2)
 279     nb_kernel_setup_ia32_sse2(fplog,nb_kernel_list);
 280 #endif
 281
 282 #if defined(GMX_X86_64_SSE)
 283     nb_kernel_setup_x86_64_sse(fplog,nb_kernel_list);
 284 #endif
 285
 286 #if defined(GMX_X86_64_SSE2)
 287     nb_kernel_setup_x86_64_sse2(fplog,nb_kernel_list);
 288 #endif
 289
 290 #if (defined GMX_IA64_ASM && defined GMX_DOUBLE)
 291     nb_kernel_setup_ia64_double(fplog,nb_kernel_list);
 292 #endif
 293
 294 #if (defined GMX_IA64_ASM && !defined GMX_DOUBLE)
 295     nb_kernel_setup_ia64_single(fplog,nb_kernel_list);
 296 #endif
 297
 298         if(fplog)
 299     {
 300             fprintf(fplog,"\n\n");
 301     }
 302 }
 303
 304 void
 305 gmx_setup_adress_kernels(FILE *fplog,gmx_bool bGenericKernelOnly) {
 306     int i;
 307
 308     snew(nb_kernel_list_adress, eNR_NBKERNEL_NR);
 309
 310     for (i = 0; i < eNR_NBKERNEL_NR; i++) {
 311         nb_kernel_list_adress[i] = NULL;
 312     }
 313
 314     if (bGenericKernelOnly)
 315     {
 316         return;
 317     }
 318
 319     nb_kernel_setup_adress(fplog, nb_kernel_list_adress);
 320 }
 321
 322 void do_nonbonded(t_commrec *cr,t_forcerec *fr,
 323                   rvec x[],rvec f[],t_mdatoms *mdatoms,t_blocka *excl,
 324                   real egnb[],real egcoul[],real egpol[],rvec box_size,
 325                   t_nrnb *nrnb,real lambda,real *dvdlambda,
 326                   int nls,int eNL,int flags)
 327 {
 328     gmx_bool            bLR,bDoForces,bForeignLambda;
 329         t_nblist *      nlist;
 330         real *          fshift;
 331         int             n,n0,n1,i,i0,i1,nrnb_ind,sz;
 332         t_nblists       *nblists;
 333         gmx_bool            bWater;
 334         nb_kernel_t *   kernelptr;
 335         nb_adress_kernel_t * adresskernelptr;
 336         FILE *          fp;
 337         int             fac=0;
 338         int             nthreads = 1;
 339         int             tabletype;
 340         int             outeriter,inneriter;
 341         real *          tabledata = NULL;
 342         gmx_gbdata_t    gbdata;
 343
 344         gmx_bool        bCG; /* for AdresS */
 345         int             k;/* for AdresS */
 346
 347     bLR            = (flags & GMX_DONB_LR);
 348     bDoForces      = (flags & GMX_DONB_FORCES);
 349     bForeignLambda = (flags & GMX_DONB_FOREIGNLAMBDA);
 350
 351     bCG = FALSE;  /* for AdresS */
 352     adresskernelptr = NULL;
 353
 354         gbdata.gb_epsilon_solvent = fr->gb_epsilon_solvent;
 355         gbdata.epsilon_r = fr->epsilon_r;
 356         gbdata.gpol               = egpol;
 357
 358     if (!fr->adress_type==eAdressOff && !bDoForces){
 359         gmx_fatal(FARGS,"No force kernels not implemeted for adress");
 360     }
 361
 362     if(fr->bAllvsAll)
 363     {
 364         if(fr->bGB)
 365         {
 366 #if (defined GMX_SSE2 || defined GMX_X86_64_SSE || defined GMX_X86_64_SSE2 || defined GMX_IA32_SSE || defined GMX_IA32_SSE2)
 367 # ifdef GMX_DOUBLE
 368             if(fr->UseOptimizedKernels)
 369             {
 370                 nb_kernel_allvsallgb_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
 371                                                  &outeriter,&inneriter,&fr->AllvsAll_work);
 372             }
 373             else
 374             {
 375                 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
 376                                      &outeriter,&inneriter,&fr->AllvsAll_work);
 377             }
 378 #  else /* not double */
 379             if(fr->UseOptimizedKernels)
 380             {
 381                 nb_kernel_allvsallgb_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
 382                                                  &outeriter,&inneriter,&fr->AllvsAll_work);
 383             }
 384             else
 385             {
 386                 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
 387                                      &outeriter,&inneriter,&fr->AllvsAll_work);
 388             }
 389 #  endif /* double/single alt. */
 390 #else /* no SSE support compiled in */
 391             nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
 392                                  &outeriter,&inneriter,&fr->AllvsAll_work);
 393 #endif
 394             inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALLGB,inneriter);
 395         }
 396         else
 397         {
 398 #if (defined GMX_SSE2 || defined GMX_X86_64_SSE || defined GMX_X86_64_SSE2 || defined GMX_IA32_SSE || defined GMX_IA32_SSE2)
 399 # ifdef GMX_DOUBLE
 400             if(fr->UseOptimizedKernels)
 401             {
 402                 nb_kernel_allvsall_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
 403                                                &outeriter,&inneriter,&fr->AllvsAll_work);
 404             }
 405             else
 406             {
 407                 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
 408                                    &outeriter,&inneriter,&fr->AllvsAll_work);
 409             }
 410
 411 #  else /* not double */
 412             if(fr->UseOptimizedKernels)
 413             {
 414                 nb_kernel_allvsall_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
 415                                                &outeriter,&inneriter,&fr->AllvsAll_work);
 416             }
 417             else
 418             {
 419                 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
 420                                    &outeriter,&inneriter,&fr->AllvsAll_work);
 421             }
 422
 423 #  endif /* double/single check */
 424 #else /* No SSE2 support compiled in */
 425             nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
 426                                &outeriter,&inneriter,&fr->AllvsAll_work);
 427 #endif
 428
 429             inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALL,inneriter);
 430         }
 431         inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,outeriter);
 432         return;
 433     }
 434
 435     if (eNL >= 0)
 436     {
 437                 i0 = eNL;
 438                 i1 = i0+1;
 439     }
 440     else
 441     {
 442                 i0 = 0;
 443                 i1 = eNL_NR;
 444         }
 445
 446         if (nls >= 0)
 447         {
 448                 n0 = nls;
 449                 n1 = nls+1;
 450         }
 451         else
 452         {
 453                 n0 = 0;
 454                 n1 = fr->nnblists;
 455         }
 456
 457         if(nb_kernel_list == NULL)
 458     {
 459                 gmx_fatal(FARGS,"gmx_setup_kernels has not been called");
 460     }
 461
 462     fshift = fr->fshift[0];
 463
 464         for(n=n0; (n<n1); n++)
 465         {
 466                 nblists = &fr->nblists[n];
 467                 for(i=i0; (i<i1); i++)
 468                 {
 469                         outeriter = inneriter = 0;
 470
 471                         if (bLR)
 472                         {
 473                                 nlist = &(nblists->nlist_lr[i]);
 474                         }
 475                         else
 476                         {
 477                                 nlist = &(nblists->nlist_sr[i]);
 478                         }
 479
 480                         if (nlist->nri > 0)
 481                         {
 482                                 nrnb_ind = nlist->il_code;
 483
 484                                 if(nrnb_ind==eNR_NBKERNEL_FREE_ENERGY)
 485                                 {
 486                                         /* generic free energy, use combined table */
 487                                         tabledata = nblists->tab.tab;
 488                                 }
 489                                 else
 490                                 {
 491                     if (bForeignLambda)
 492                     {
 493                         /* We don't need the non-perturbed interactions */
 494                         continue;
 495                     }
 496
 497                                         tabletype = nb_kernel_table[nrnb_ind];
 498
 499                                         /* normal kernels, not free energy */
 500                                         if (!bDoForces)
 501                                         {
 502                                                 nrnb_ind += eNR_NBKERNEL_NR/2;
 503                                         }
 504
 505                                         if(tabletype == TABLE_COMBINED)
 506                                         {
 507                                                 tabledata = nblists->tab.tab;
 508                                         }
 509                                         else if(tabletype == TABLE_COUL)
 510                                         {
 511                                                 tabledata = nblists->coultab;
 512                                         }
 513                                         else if(tabletype == TABLE_VDW)
 514                                         {
 515                                                 tabledata = nblists->vdwtab;
 516                                         }
 517                                         else
 518                                         {
 519                                                 tabledata = NULL;
 520                                         }
 521                                 }
 522
 523                                 nlist->count = 0;
 524
 525                                 if(nlist->free_energy)
 526                                 {
 527                                         if(nlist->ivdw==2)
 528                                         {
 529                                                 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
 530                                         }
 531
 532                                         gmx_nb_free_energy_kernel(nlist->icoul,
 533                                                                                           nlist->ivdw,
 534                                                                                           nlist->nri,
 535                                                                                           nlist->iinr,
 536                                                                                           nlist->jindex,
 537                                                                                           nlist->jjnr,
 538                                                                                           nlist->shift,
 539                                                                                           fr->shift_vec[0],
 540                                                                                           fshift,
 541                                                                                           nlist->gid,
 542                                                                                           x[0],
 543                                                                                           f[0],
 544                                                                                           mdatoms->chargeA,
 545                                                                                           mdatoms->chargeB,
 546                                                                                           fr->epsfac,
 547                                                                                           fr->k_rf,
 548                                                                                           fr->c_rf,
 549                                                                                           fr->ewaldcoeff,
 550                                                                                           egcoul,
 551                                                                                           mdatoms->typeA,
 552                                                                                           mdatoms->typeB,
 553                                                                                           fr->ntype,
 554                                                                                           fr->nbfp,
 555                                                                                           egnb,
 556                                                                                           nblists->tab.scale,
 557                                                                                           tabledata,
 558                                                                                           lambda,
 559                                                                                           dvdlambda,
 560                                                                                           fr->sc_alpha,
 561                                                                                           fr->sc_power,
 562                                                                                           fr->sc_sigma6_def,
 563                                               fr->sc_sigma6_min,
 564                                               bDoForces,
 565                                                                                           &outeriter,
 566                                                                                           &inneriter);
 567                 }
 568                 else if (nlist->enlist == enlistCG_CG)
 569                 {
 570                     if (fr->adress_type==eAdressOff){
 571                     /* Call the charge group based inner loop */
 572                        gmx_nb_generic_cg_kernel(nlist,
 573                                                 fr,
 574                                                 mdatoms,
 575                                                 x[0],
 576                                                 f[0],
 577                                                 fshift,
 578                                                 egcoul,
 579                                                 egnb,
 580                                                 nblists->tab.scale,
 581                                                 tabledata,
 582                                                 &outeriter,
 583                                                 &inneriter);
 584                     }
 585                     else
 586                     {
 587                        /*gmx_nb_generic_adress_kernel(nlist,
 588                                                 fr,
 589                                                 mdatoms,
 590                                                 x[0],
 591                                                 f[0],
 592                                                 fshift,
 593                                                 egcoul,
 594                                                 egnb,
 595                                                 nblists->tab.scale,
 596                                                 tabledata,
 597                                                 &outeriter,
 598                                                 &inneriter);*/
 599                           gmx_fatal(FARGS,"Death & horror! Adress cgcg kernel not implemented anymore.\n");
 600
 601                     }
 602                 }
 603                 else
 604                 {
 605                     /* AdresS*/
 606                     /* for adress we need to determine for each energy group wether it is explicit or coarse-grained */
 607                     if (!fr->adress_type == eAdressOff) {
 608                         bCG = FALSE;
 609                         if ( !fr->adress_group_explicit[ mdatoms->cENER[nlist->iinr[0]] ] ){
 610                             bCG=TRUE;
 611                         }
 612                         /* If this processor has only explicit atoms (w=1)
 613                           skip the coarse grained force calculation. Same for
 614                          only coarsegrained atoms and explicit interactions.
 615                          Last condition is to make sure that generic kernel is not
 616                          skipped*/
 617                         if (mdatoms->pureex && bCG && nb_kernel_list[nrnb_ind] != NULL) continue;
 618                         if (mdatoms->purecg && !bCG && nb_kernel_list[nrnb_ind] != NULL) continue;
 619                         kernelptr = NULL;
 620                         adresskernelptr = NULL;
 621                     }
 622
 623                     if (fr->adress_type == eAdressOff ||
 624                             mdatoms->pureex ||
 625                             mdatoms->purecg){
 626                         /* if we only have to calculate pure cg/ex interactions
 627                          we can use the faster standard gromacs kernels*/
 628                         kernelptr = nb_kernel_list[nrnb_ind];
 629                     }else{
 630                         /* This processor has hybrid interactions which means
 631                          * we have to
 632                          * use our own kernels. We have two kernel types: one that
 633                          * calculates the forces with the explicit prefactor w1*w2
 634                          * and one for coarse-grained with (1-w1*w2)
 635                          * explicit kernels are the second part of the kernel
 636                          *  list */
 637                         if (!bCG) nrnb_ind += eNR_NBKERNEL_NR/2;
 638                         adresskernelptr = nb_kernel_list_adress[nrnb_ind];
 639                     }
 640
 641                     if (kernelptr == NULL && adresskernelptr == NULL)
 642                      {
 643                         /* Call a generic nonbonded kernel */
 644
 645                         /* If you want to hack/test your own interactions,
 646                          * do it in this routine and make sure it is called
 647                          * by setting the environment variable GMX_NB_GENERIC.
 648                          */
 649                         if (fr->adress_type==eAdressOff){
 650
 651                         gmx_nb_generic_kernel(nlist,
 652                                               fr,
 653                                               mdatoms,
 654                                               x[0],
 655                                               f[0],
 656                                               fshift,
 657                                               egcoul,
 658                                               egnb,
 659                                               nblists->tab.scale,
 660                                               tabledata,
 661                                               &outeriter,
 662                                               &inneriter);
 663                         }else /* do generic AdResS kernels (slow)*/
 664                         {
 665
 666                             gmx_nb_generic_adress_kernel(nlist,
 667                                                 fr,
 668                                                 mdatoms,
 669                                                 x[0],
 670                                                 f[0],
 671                                                 fshift,
 672                                                 egcoul,
 673                                                 egnb,
 674                                                 nblists->tab.scale,
 675                                                 tabledata,
 676                                                 &outeriter,
 677                                                 &inneriter,
 678                                                 bCG);
 679                         }
 680
 681
 682                     }
 683                     else
 684                     {
 685                         /* Call nonbonded kernel from function pointer */
 686                         if (kernelptr!=NULL){
 687                         (*kernelptr)( &(nlist->nri),
 688                                       nlist->iinr,
 689                                       nlist->jindex,
 690                                       nlist->jjnr,
 691                                       nlist->shift,
 692                                       fr->shift_vec[0],
 693                                       fshift,
 694                                       nlist->gid,
 695                                       x[0],
 696                                       f[0],
 697                                       mdatoms->chargeA,
 698                                       &(fr->epsfac),
 699                                       &(fr->k_rf),
 700                                       &(fr->c_rf),
 701                                       egcoul,
 702                                       mdatoms->typeA,
 703                                       &(fr->ntype),
 704                                       fr->nbfp,
 705                                       egnb,
 706                                       &(nblists->tab.scale),
 707                                       tabledata,
 708                                       fr->invsqrta,
 709                                       fr->dvda,
 710                                       &(fr->gbtabscale),
 711                                       fr->gbtab.tab,
 712                                       &nthreads,
 713                                       &(nlist->count),
 714                                       nlist->mtx,
 715                                       &outeriter,
 716                                       &inneriter,
 717                                       (real *)&gbdata);
 718                         }else if (adresskernelptr != NULL)
 719                         { /* Adress kernels */
 720                           (*adresskernelptr)( &(nlist->nri),
 721                                       nlist->iinr,
 722                                       nlist->jindex,
 723                                       nlist->jjnr,
 724                                       nlist->shift,
 725                                       fr->shift_vec[0],
 726                                       fshift,
 727                                       nlist->gid,
 728                                       x[0],
 729                                       f[0],
 730                                       mdatoms->chargeA,
 731                                       &(fr->epsfac),
 732                                       &(fr->k_rf),
 733                                       &(fr->c_rf),
 734                                       egcoul,
 735                                       mdatoms->typeA,
 736                                       &(fr->ntype),
 737                                       fr->nbfp,
 738                                       egnb,
 739                                       &(nblists->tab.scale),
 740                                       tabledata,
 741                                       fr->invsqrta,
 742                                       fr->dvda,
 743                                       &(fr->gbtabscale),
 744                                       fr->gbtab.tab,
 745                                       &nthreads,
 746                                       &(nlist->count),
 747                                       nlist->mtx,
 748                                       &outeriter,
 749                                       &inneriter,
 750                                       fr->adress_ex_forcecap,
 751                                       mdatoms->wf);
 752                         }
 753                     }
 754                 }
 755
 756                 /* Update flop accounting */
 757
 758                                 /* Outer loop in kernel */
 759                 switch (nlist->enlist) {
 760                 case enlistATOM_ATOM:   fac =  1; break;
 761                 case enlistSPC_ATOM:    fac =  3; break;
 762                 case enlistSPC_SPC:     fac =  9; break;
 763                 case enlistTIP4P_ATOM:  fac =  4; break;
 764                 case enlistTIP4P_TIP4P: fac = 16; break;
 765                 case enlistCG_CG:       fac =  1; break;
 766                 }
 767                 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fac*outeriter);
 768
 769                 /* inner loop in kernel */
 770                 inc_nrnb(nrnb,nrnb_ind,inneriter);
 771             }
 772         }
 773     }
 774 }
 775
 776
 777 real
 778 do_listed_vdw_q(int ftype,int nbonds,
 779                 const t_iatom iatoms[],const t_iparams iparams[],
 780                 const rvec x[],rvec f[],rvec fshift[],
 781                 const t_pbc *pbc,const t_graph *g,
 782                 real lambda,real *dvdlambda,
 783                 const t_mdatoms *md,
 784                 const t_forcerec *fr,gmx_grppairener_t *grppener,
 785                 int *global_atom_index)
 786 {
 787     static    gmx_bool bWarn=FALSE;
 788     real      eps,r2,*tab,rtab2=0;
 789     rvec      dx,x14[2],f14[2];
 790     int       i,ai,aj,itype;
 791     int       typeA[2]={0,0},typeB[2]={0,1};
 792     real      chargeA[2]={0,0},chargeB[2];
 793     int       gid,shift_vir,shift_f;
 794     int       j_index[] = { 0, 1 };
 795     int       i0=0,i1=1,i2=2;
 796     ivec      dt;
 797     int       outeriter,inneriter;
 798     int       nthreads = 1;
 799     int       count;
 800     real      krf,crf,tabscale;
 801     int       ntype=0;
 802     real      *nbfp=NULL;
 803     real      *egnb=NULL,*egcoul=NULL;
 804     t_nblist  tmplist;
 805     int       icoul,ivdw;
 806     gmx_bool      bMolPBC,bFreeEnergy;
 807
 808     gmx_bool      bCG; /* AdResS*/
 809     real      wf14[2]={0,0}; /* AdResS*/
 810
 811 #if GMX_THREAD_SHM_FDECOMP
 812     pthread_mutex_t mtx;
 813 #else
 814     void *    mtx = NULL;
 815 #endif
 816
 817
 818 #if GMX_THREAD_SHM_FDECOMP
 819     pthread_mutex_initialize(&mtx);
 820 #endif
 821
 822     bMolPBC = fr->bMolPBC;
 823
 824     switch (ftype) {
 825     case F_LJ14:
 826     case F_LJC14_Q:
 827         eps = fr->epsfac*fr->fudgeQQ;
 828         ntype  = 1;
 829         egnb   = grppener->ener[egLJ14];
 830         egcoul = grppener->ener[egCOUL14];
 831         break;
 832     case F_LJC_PAIRS_NB:
 833         eps = fr->epsfac;
 834         ntype  = 1;
 835         egnb   = grppener->ener[egLJSR];
 836         egcoul = grppener->ener[egCOULSR];
 837         break;
 838     default:
 839         gmx_fatal(FARGS,"Unknown function type %d in do_nonbonded14",
 840                   ftype);
 841     }
 842     tab = fr->tab14.tab;
 843     rtab2 = sqr(fr->tab14.r);
 844     tabscale = fr->tab14.scale;
 845
 846     krf = fr->k_rf;
 847     crf = fr->c_rf;
 848
 849     /* Determine the values for icoul/ivdw. */
 850     if (fr->bEwald) {
 851         icoul = 1;
 852     }
 853     else if(fr->bcoultab)
 854     {
 855         icoul = 3;
 856     }
 857     else if(fr->eeltype == eelRF_NEC)
 858     {
 859         icoul = 2;
 860     }
 861     else
 862     {
 863         icoul = 1;
 864     }
 865
 866     if(fr->bvdwtab)
 867     {
 868         ivdw = 3;
 869     }
 870     else if(fr->bBHAM)
 871     {
 872         ivdw = 2;
 873     }
 874     else
 875     {
 876         ivdw = 1;
 877     }
 878
 879
 880     bCG = FALSE; /*Adres*/
 881     /* We don't do SSE or altivec here, due to large overhead for 4-fold
 882      * unrolling on short lists
 883      */
 884
 885     bFreeEnergy = FALSE;
 886     for(i=0; (i<nbonds); )
 887     {
 888         itype = iatoms[i++];
 889         ai    = iatoms[i++];
 890         aj    = iatoms[i++];
 891         gid   = GID(md->cENER[ai],md->cENER[aj],md->nenergrp);
 892
 893         if (!fr->adress_type == eAdressOff) {
 894             if (fr->adress_group_explicit[md->cENER[ai]] != fr->adress_group_explicit[md->cENER[aj]]){
 895                 /*exclude cg-ex interaction*/
 896                 continue;
 897             }
 898             bCG = !fr->adress_group_explicit[md->cENER[ai]];
 899             wf14[0] = md->wf[ai];
 900             wf14[1] = md->wf[aj];
 901         }
 902         switch (ftype) {
 903         case F_LJ14:
 904             bFreeEnergy =
 905                 (fr->efep != efepNO &&
 906                  ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
 907                   iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
 908                   iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 909             chargeA[0] = md->chargeA[ai];
 910             chargeA[1] = md->chargeA[aj];
 911             nbfp = (real *)&(iparams[itype].lj14.c6A);
 912             break;
 913         case F_LJC14_Q:
 914             eps = fr->epsfac*iparams[itype].ljc14.fqq;
 915             chargeA[0] = iparams[itype].ljc14.qi;
 916             chargeA[1] = iparams[itype].ljc14.qj;
 917             nbfp = (real *)&(iparams[itype].ljc14.c6);
 918             break;
 919         case F_LJC_PAIRS_NB:
 920             chargeA[0] = iparams[itype].ljcnb.qi;
 921             chargeA[1] = iparams[itype].ljcnb.qj;
 922             nbfp = (real *)&(iparams[itype].ljcnb.c6);
 923             break;
 924         }
 925
 926         if (!bMolPBC)
 927         {
 928             /* This is a bonded interaction, atoms are in the same box */
 929             shift_f = CENTRAL;
 930             r2 = distance2(x[ai],x[aj]);
 931         }
 932         else
 933         {
 934             /* Apply full periodic boundary conditions */
 935             shift_f = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
 936             r2 = norm2(dx);
 937         }
 938
 939         if (r2 >= rtab2)
 940         {
 941             if (!bWarn)
 942             {
 943                 fprintf(stderr,"Warning: 1-4 interaction between %d and %d "
 944                         "at distance %.3f which is larger than the 1-4 table size %.3f nm\n",
 945                         glatnr(global_atom_index,ai),
 946                         glatnr(global_atom_index,aj),
 947                         sqrt(r2), sqrt(rtab2));
 948                 fprintf(stderr,"These are ignored for the rest of the simulation\n");
 949                 fprintf(stderr,"This usually means your system is exploding,\n"
 950                         "if not, you should increase table-extension in your mdp file\n"
 951                         "or with user tables increase the table size\n");
 952                 bWarn = TRUE;
 953             }
 954             if (debug)
 955               fprintf(debug,"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
 956                       x[ai][XX],x[ai][YY],x[ai][ZZ],
 957                       x[aj][XX],x[aj][YY],x[aj][ZZ],
 958                       glatnr(global_atom_index,ai),
 959                       glatnr(global_atom_index,aj),
 960                       sqrt(r2));
 961         }
 962         else
 963         {
 964             copy_rvec(x[ai],x14[0]);
 965             copy_rvec(x[aj],x14[1]);
 966             clear_rvec(f14[0]);
 967             clear_rvec(f14[1]);
 968 #ifdef DEBUG
 969             fprintf(debug,"LJ14: grp-i=%2d, grp-j=%2d, ngrp=%2d, GID=%d\n",
 970                     md->cENER[ai],md->cENER[aj],md->nenergrp,gid);
 971 #endif
 972
 973             outeriter = inneriter = count = 0;
 974             if (bFreeEnergy)
 975         {
 976             chargeB[0] = md->chargeB[ai];
 977             chargeB[1] = md->chargeB[aj];
 978             /* We pass &(iparams[itype].lj14.c6A) as LJ parameter matrix
 979              * to the innerloops.
 980              * Here we use that the LJ-14 parameters are stored in iparams
 981              * as c6A,c12A,c6B,c12B, which are referenced correctly
 982              * in the innerloops if we assign type combinations 0-0 and 0-1
 983              * to atom pair ai-aj in topologies A and B respectively.
 984              */
 985             if(ivdw==2)
 986             {
 987                 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
 988             }
 989             count = 0;
 990             gmx_nb_free_energy_kernel(icoul,
 991                                       ivdw,
 992                                       i1,
 993                                       &i0,
 994                                       j_index,
 995                                       &i1,
 996                                       &shift_f,
 997                                       fr->shift_vec[0],
 998                                       fshift[0],
 999                                       &gid,
1000                                       x14[0],
1001                                       f14[0],
1002                                       chargeA,
1003                                       chargeB,
1004                                       eps,
1005                                       krf,
1006                                       crf,
1007                                       fr->ewaldcoeff,
1008                                       egcoul,
1009                                       typeA,
1010                                       typeB,
1011                                       ntype,
1012                                       nbfp,
1013                                       egnb,
1014                                       tabscale,
1015                                       tab,
1016                                       lambda,
1017                                       dvdlambda,
1018                                       fr->sc_alpha,
1019                                       fr->sc_power,
1020                                       fr->sc_sigma6_def,
1021                                       fr->sc_sigma6_min,
1022                                       TRUE,
1023                                       &outeriter,
1024                                       &inneriter);
1025         }
1026         else
1027         {
1028           if (fr->adress_type==eAdressOff || !fr->adress_do_hybridpairs){
1029             /* Not perturbed - call kernel 330 */
1030             nb_kernel330
1031                 ( &i1,
1032                   &i0,
1033                   j_index,
1034                   &i1,
1035                   &shift_f,
1036                   fr->shift_vec[0],
1037                   fshift[0],
1038                   &gid,
1039                   x14[0],
1040                   f14[0],
1041                   chargeA,
1042                   &eps,
1043                   &krf,
1044                   &crf,
1045                   egcoul,
1046                   typeA,
1047                   &ntype,
1048                   nbfp,
1049                   egnb,
1050                   &tabscale,
1051                   tab,
1052                   NULL,
1053                   NULL,
1054                   NULL,
1055                   NULL,
1056                   &nthreads,
1057                   &count,
1058                   (void *)&mtx,
1059                   &outeriter,
1060                   &inneriter,
1061                   NULL);
1062                 } else {
1063                     if (bCG) {
1064                         nb_kernel330_adress_cg(&i1,
1065                                 &i0,
1066                                 j_index,
1067                                 &i1,
1068                                 &shift_f,
1069                                 fr->shift_vec[0],
1070                                 fshift[0],
1071                                 &gid,
1072                                 x14[0],
1073                                 f14[0],
1074                                 chargeA,
1075                                 &eps,
1076                                 &krf,
1077                                 &crf,
1078                                 egcoul,
1079                                 typeA,
1080                                 &ntype,
1081                                 nbfp,
1082                                 egnb,
1083                                 &tabscale,
1084                                 tab,
1085                                 NULL,
1086                                 NULL,
1087                                 NULL,
1088                                 NULL,
1089                                 &nthreads,
1090                                 &count,
1091                                 (void *) &mtx,
1092                                 &outeriter,
1093                                 &inneriter,
1094                                 fr->adress_ex_forcecap,
1095                                 wf14);
1096                     } else {
1097                         nb_kernel330_adress_ex(&i1,
1098                                 &i0,
1099                                 j_index,
1100                                 &i1,
1101                                 &shift_f,
1102                                 fr->shift_vec[0],
1103                                 fshift[0],
1104                                 &gid,
1105                                 x14[0],
1106                                 f14[0],
1107                                 chargeA,
1108                                 &eps,
1109                                 &krf,
1110                                 &crf,
1111                                 egcoul,
1112                                 typeA,
1113                                 &ntype,
1114                                 nbfp,
1115                                 egnb,
1116                                 &tabscale,
1117                                 tab,
1118                                 NULL,
1119                                 NULL,
1120                                 NULL,
1121                                 NULL,
1122                                 &nthreads,
1123                                 &count,
1124                                 (void *) &mtx,
1125                                 &outeriter,
1126                                 &inneriter,
1127                                 fr->adress_ex_forcecap,
1128                                 wf14);
1129                     }
1130
1131                 }
1132             }
1133
1134         /* Add the forces */
1135         rvec_inc(f[ai],f14[0]);
1136         rvec_dec(f[aj],f14[0]);
1137
1138         if (g)
1139         {
1140             /* Correct the shift forces using the graph */
1141             ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1142             shift_vir = IVEC2IS(dt);
1143             rvec_inc(fshift[shift_vir],f14[0]);
1144             rvec_dec(fshift[CENTRAL],f14[0]);
1145         }
1146
1147             /* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
1148         }
1149     }
1150     return 0.0;
1151 }
1152
1153