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         FILE *          fp;
 335         int             fac=0;
 336         int             nthreads = 1;
 337         int             tabletype;
 338         int             outeriter,inneriter;
 339         real *          tabledata = NULL;
 340         gmx_gbdata_t    gbdata;
 341
 342         nb_kernel_t         *kernelptr=NULL;
 343         nb_adress_kernel_t  *adresskernelptr=NULL;
 344
 345         gmx_bool        bCG; /* for AdresS */
 346         int             k;/* for AdresS */
 347
 348     bLR            = (flags & GMX_DONB_LR);
 349     bDoForces      = (flags & GMX_DONB_FORCES);
 350     bForeignLambda = (flags & GMX_DONB_FOREIGNLAMBDA);
 351
 352     bCG = FALSE;  /* for AdresS */
 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                     }
 620
 621                     if (fr->adress_type == eAdressOff ||
 622                             mdatoms->pureex ||
 623                             mdatoms->purecg){
 624                         /* if we only have to calculate pure cg/ex interactions
 625                          we can use the faster standard gromacs kernels*/
 626                         kernelptr = nb_kernel_list[nrnb_ind];
 627                     }else{
 628                         /* This processor has hybrid interactions which means
 629                          * we have to
 630                          * use our own kernels. We have two kernel types: one that
 631                          * calculates the forces with the explicit prefactor w1*w2
 632                          * and one for coarse-grained with (1-w1*w2)
 633                          * explicit kernels are the second part of the kernel
 634                          *  list */
 635                         if (!bCG) nrnb_ind += eNR_NBKERNEL_NR/2;
 636                         adresskernelptr = nb_kernel_list_adress[nrnb_ind];
 637                     }
 638
 639                     if (kernelptr == NULL && adresskernelptr == NULL)
 640                      {
 641                         /* Call a generic nonbonded kernel */
 642
 643                         /* If you want to hack/test your own interactions,
 644                          * do it in this routine and make sure it is called
 645                          * by setting the environment variable GMX_NB_GENERIC.
 646                          */
 647                         if (fr->adress_type==eAdressOff){
 648
 649                         gmx_nb_generic_kernel(nlist,
 650                                               fr,
 651                                               mdatoms,
 652                                               x[0],
 653                                               f[0],
 654                                               fshift,
 655                                               egcoul,
 656                                               egnb,
 657                                               nblists->tab.scale,
 658                                               tabledata,
 659                                               &outeriter,
 660                                               &inneriter);
 661                         }else /* do generic AdResS kernels (slow)*/
 662                         {
 663
 664                             gmx_nb_generic_adress_kernel(nlist,
 665                                                 fr,
 666                                                 mdatoms,
 667                                                 x[0],
 668                                                 f[0],
 669                                                 fshift,
 670                                                 egcoul,
 671                                                 egnb,
 672                                                 nblists->tab.scale,
 673                                                 tabledata,
 674                                                 &outeriter,
 675                                                 &inneriter,
 676                                                 bCG);
 677                         }
 678
 679
 680                     }
 681                     else
 682                     {
 683                         /* Call nonbonded kernel from function pointer */
 684                         if (kernelptr!=NULL){
 685                         (*kernelptr)( &(nlist->nri),
 686                                       nlist->iinr,
 687                                       nlist->jindex,
 688                                       nlist->jjnr,
 689                                       nlist->shift,
 690                                       fr->shift_vec[0],
 691                                       fshift,
 692                                       nlist->gid,
 693                                       x[0],
 694                                       f[0],
 695                                       mdatoms->chargeA,
 696                                       &(fr->epsfac),
 697                                       &(fr->k_rf),
 698                                       &(fr->c_rf),
 699                                       egcoul,
 700                                       mdatoms->typeA,
 701                                       &(fr->ntype),
 702                                       fr->nbfp,
 703                                       egnb,
 704                                       &(nblists->tab.scale),
 705                                       tabledata,
 706                                       fr->invsqrta,
 707                                       fr->dvda,
 708                                       &(fr->gbtabscale),
 709                                       fr->gbtab.tab,
 710                                       &nthreads,
 711                                       &(nlist->count),
 712                                       nlist->mtx,
 713                                       &outeriter,
 714                                       &inneriter,
 715                                       (real *)&gbdata);
 716                         }else if (adresskernelptr != NULL)
 717                         { /* Adress kernels */
 718                           (*adresskernelptr)( &(nlist->nri),
 719                                       nlist->iinr,
 720                                       nlist->jindex,
 721                                       nlist->jjnr,
 722                                       nlist->shift,
 723                                       fr->shift_vec[0],
 724                                       fshift,
 725                                       nlist->gid,
 726                                       x[0],
 727                                       f[0],
 728                                       mdatoms->chargeA,
 729                                       &(fr->epsfac),
 730                                       &(fr->k_rf),
 731                                       &(fr->c_rf),
 732                                       egcoul,
 733                                       mdatoms->typeA,
 734                                       &(fr->ntype),
 735                                       fr->nbfp,
 736                                       egnb,
 737                                       &(nblists->tab.scale),
 738                                       tabledata,
 739                                       fr->invsqrta,
 740                                       fr->dvda,
 741                                       &(fr->gbtabscale),
 742                                       fr->gbtab.tab,
 743                                       &nthreads,
 744                                       &(nlist->count),
 745                                       nlist->mtx,
 746                                       &outeriter,
 747                                       &inneriter,
 748                                       fr->adress_ex_forcecap,
 749                                       mdatoms->wf);
 750                         }
 751                     }
 752                 }
 753
 754                 /* Update flop accounting */
 755
 756                                 /* Outer loop in kernel */
 757                 switch (nlist->enlist) {
 758                 case enlistATOM_ATOM:   fac =  1; break;
 759                 case enlistSPC_ATOM:    fac =  3; break;
 760                 case enlistSPC_SPC:     fac =  9; break;
 761                 case enlistTIP4P_ATOM:  fac =  4; break;
 762                 case enlistTIP4P_TIP4P: fac = 16; break;
 763                 case enlistCG_CG:       fac =  1; break;
 764                 }
 765                 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fac*outeriter);
 766
 767                 /* inner loop in kernel */
 768                 inc_nrnb(nrnb,nrnb_ind,inneriter);
 769             }
 770         }
 771     }
 772 }
 773
 774
 775 real
 776 do_listed_vdw_q(int ftype,int nbonds,
 777                 const t_iatom iatoms[],const t_iparams iparams[],
 778                 const rvec x[],rvec f[],rvec fshift[],
 779                 const t_pbc *pbc,const t_graph *g,
 780                 real lambda,real *dvdlambda,
 781                 const t_mdatoms *md,
 782                 const t_forcerec *fr,gmx_grppairener_t *grppener,
 783                 int *global_atom_index)
 784 {
 785     static    gmx_bool bWarn=FALSE;
 786     real      eps,r2,*tab,rtab2=0;
 787     rvec      dx,x14[2],f14[2];
 788     int       i,ai,aj,itype;
 789     int       typeA[2]={0,0},typeB[2]={0,1};
 790     real      chargeA[2]={0,0},chargeB[2];
 791     int       gid,shift_vir,shift_f;
 792     int       j_index[] = { 0, 1 };
 793     int       i0=0,i1=1,i2=2;
 794     ivec      dt;
 795     int       outeriter,inneriter;
 796     int       nthreads = 1;
 797     int       count;
 798     real      krf,crf,tabscale;
 799     int       ntype=0;
 800     real      *nbfp=NULL;
 801     real      *egnb=NULL,*egcoul=NULL;
 802     t_nblist  tmplist;
 803     int       icoul,ivdw;
 804     gmx_bool      bMolPBC,bFreeEnergy;
 805
 806     gmx_bool      bCG; /* AdResS*/
 807     real      wf14[2]={0,0}; /* AdResS*/
 808
 809 #if GMX_THREAD_SHM_FDECOMP
 810     pthread_mutex_t mtx;
 811 #else
 812     void *    mtx = NULL;
 813 #endif
 814
 815
 816 #if GMX_THREAD_SHM_FDECOMP
 817     pthread_mutex_initialize(&mtx);
 818 #endif
 819
 820     bMolPBC = fr->bMolPBC;
 821
 822     switch (ftype) {
 823     case F_LJ14:
 824     case F_LJC14_Q:
 825         eps = fr->epsfac*fr->fudgeQQ;
 826         ntype  = 1;
 827         egnb   = grppener->ener[egLJ14];
 828         egcoul = grppener->ener[egCOUL14];
 829         break;
 830     case F_LJC_PAIRS_NB:
 831         eps = fr->epsfac;
 832         ntype  = 1;
 833         egnb   = grppener->ener[egLJSR];
 834         egcoul = grppener->ener[egCOULSR];
 835         break;
 836     default:
 837         gmx_fatal(FARGS,"Unknown function type %d in do_nonbonded14",
 838                   ftype);
 839     }
 840     tab = fr->tab14.tab;
 841     rtab2 = sqr(fr->tab14.r);
 842     tabscale = fr->tab14.scale;
 843
 844     krf = fr->k_rf;
 845     crf = fr->c_rf;
 846
 847     /* Determine the values for icoul/ivdw. */
 848     if (fr->bEwald) {
 849         icoul = 1;
 850     }
 851     else if(fr->bcoultab)
 852     {
 853         icoul = 3;
 854     }
 855     else if(fr->eeltype == eelRF_NEC)
 856     {
 857         icoul = 2;
 858     }
 859     else
 860     {
 861         icoul = 1;
 862     }
 863
 864     if(fr->bvdwtab)
 865     {
 866         ivdw = 3;
 867     }
 868     else if(fr->bBHAM)
 869     {
 870         ivdw = 2;
 871     }
 872     else
 873     {
 874         ivdw = 1;
 875     }
 876
 877
 878     bCG = FALSE; /*Adres*/
 879     /* We don't do SSE or altivec here, due to large overhead for 4-fold
 880      * unrolling on short lists
 881      */
 882
 883     bFreeEnergy = FALSE;
 884     for(i=0; (i<nbonds); )
 885     {
 886         itype = iatoms[i++];
 887         ai    = iatoms[i++];
 888         aj    = iatoms[i++];
 889         gid   = GID(md->cENER[ai],md->cENER[aj],md->nenergrp);
 890
 891         if (!fr->adress_type == eAdressOff) {
 892             if (fr->adress_group_explicit[md->cENER[ai]] != fr->adress_group_explicit[md->cENER[aj]]){
 893                 /*exclude cg-ex interaction*/
 894                 continue;
 895             }
 896             bCG = !fr->adress_group_explicit[md->cENER[ai]];
 897             wf14[0] = md->wf[ai];
 898             wf14[1] = md->wf[aj];
 899         }
 900         switch (ftype) {
 901         case F_LJ14:
 902             bFreeEnergy =
 903                 (fr->efep != efepNO &&
 904                  ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
 905                   iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
 906                   iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 907             chargeA[0] = md->chargeA[ai];
 908             chargeA[1] = md->chargeA[aj];
 909             nbfp = (real *)&(iparams[itype].lj14.c6A);
 910             break;
 911         case F_LJC14_Q:
 912             eps = fr->epsfac*iparams[itype].ljc14.fqq;
 913             chargeA[0] = iparams[itype].ljc14.qi;
 914             chargeA[1] = iparams[itype].ljc14.qj;
 915             nbfp = (real *)&(iparams[itype].ljc14.c6);
 916             break;
 917         case F_LJC_PAIRS_NB:
 918             chargeA[0] = iparams[itype].ljcnb.qi;
 919             chargeA[1] = iparams[itype].ljcnb.qj;
 920             nbfp = (real *)&(iparams[itype].ljcnb.c6);
 921             break;
 922         }
 923
 924         if (!bMolPBC)
 925         {
 926             /* This is a bonded interaction, atoms are in the same box */
 927             shift_f = CENTRAL;
 928             r2 = distance2(x[ai],x[aj]);
 929         }
 930         else
 931         {
 932             /* Apply full periodic boundary conditions */
 933             shift_f = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
 934             r2 = norm2(dx);
 935         }
 936
 937         if (r2 >= rtab2)
 938         {
 939             if (!bWarn)
 940             {
 941                 fprintf(stderr,"Warning: 1-4 interaction between %d and %d "
 942                         "at distance %.3f which is larger than the 1-4 table size %.3f nm\n",
 943                         glatnr(global_atom_index,ai),
 944                         glatnr(global_atom_index,aj),
 945                         sqrt(r2), sqrt(rtab2));
 946                 fprintf(stderr,"These are ignored for the rest of the simulation\n");
 947                 fprintf(stderr,"This usually means your system is exploding,\n"
 948                         "if not, you should increase table-extension in your mdp file\n"
 949                         "or with user tables increase the table size\n");
 950                 bWarn = TRUE;
 951             }
 952             if (debug)
 953               fprintf(debug,"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
 954                       x[ai][XX],x[ai][YY],x[ai][ZZ],
 955                       x[aj][XX],x[aj][YY],x[aj][ZZ],
 956                       glatnr(global_atom_index,ai),
 957                       glatnr(global_atom_index,aj),
 958                       sqrt(r2));
 959         }
 960         else
 961         {
 962             copy_rvec(x[ai],x14[0]);
 963             copy_rvec(x[aj],x14[1]);
 964             clear_rvec(f14[0]);
 965             clear_rvec(f14[1]);
 966 #ifdef DEBUG
 967             fprintf(debug,"LJ14: grp-i=%2d, grp-j=%2d, ngrp=%2d, GID=%d\n",
 968                     md->cENER[ai],md->cENER[aj],md->nenergrp,gid);
 969 #endif
 970
 971             outeriter = inneriter = count = 0;
 972             if (bFreeEnergy)
 973         {
 974             chargeB[0] = md->chargeB[ai];
 975             chargeB[1] = md->chargeB[aj];
 976             /* We pass &(iparams[itype].lj14.c6A) as LJ parameter matrix
 977              * to the innerloops.
 978              * Here we use that the LJ-14 parameters are stored in iparams
 979              * as c6A,c12A,c6B,c12B, which are referenced correctly
 980              * in the innerloops if we assign type combinations 0-0 and 0-1
 981              * to atom pair ai-aj in topologies A and B respectively.
 982              */
 983             if(ivdw==2)
 984             {
 985                 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
 986             }
 987             count = 0;
 988             gmx_nb_free_energy_kernel(icoul,
 989                                       ivdw,
 990                                       i1,
 991                                       &i0,
 992                                       j_index,
 993                                       &i1,
 994                                       &shift_f,
 995                                       fr->shift_vec[0],
 996                                       fshift[0],
 997                                       &gid,
 998                                       x14[0],
 999                                       f14[0],
1000                                       chargeA,
1001                                       chargeB,
1002                                       eps,
1003                                       krf,
1004                                       crf,
1005                                       fr->ewaldcoeff,
1006                                       egcoul,
1007                                       typeA,
1008                                       typeB,
1009                                       ntype,
1010                                       nbfp,
1011                                       egnb,
1012                                       tabscale,
1013                                       tab,
1014                                       lambda,
1015                                       dvdlambda,
1016                                       fr->sc_alpha,
1017                                       fr->sc_power,
1018                                       fr->sc_sigma6_def,
1019                                       fr->sc_sigma6_min,
1020                                       TRUE,
1021                                       &outeriter,
1022                                       &inneriter);
1023         }
1024         else
1025         {
1026           if (fr->adress_type==eAdressOff || !fr->adress_do_hybridpairs){
1027             /* Not perturbed - call kernel 330 */
1028             nb_kernel330
1029                 ( &i1,
1030                   &i0,
1031                   j_index,
1032                   &i1,
1033                   &shift_f,
1034                   fr->shift_vec[0],
1035                   fshift[0],
1036                   &gid,
1037                   x14[0],
1038                   f14[0],
1039                   chargeA,
1040                   &eps,
1041                   &krf,
1042                   &crf,
1043                   egcoul,
1044                   typeA,
1045                   &ntype,
1046                   nbfp,
1047                   egnb,
1048                   &tabscale,
1049                   tab,
1050                   NULL,
1051                   NULL,
1052                   NULL,
1053                   NULL,
1054                   &nthreads,
1055                   &count,
1056                   (void *)&mtx,
1057                   &outeriter,
1058                   &inneriter,
1059                   NULL);
1060                 } else {
1061                     if (bCG) {
1062                         nb_kernel330_adress_cg(&i1,
1063                                 &i0,
1064                                 j_index,
1065                                 &i1,
1066                                 &shift_f,
1067                                 fr->shift_vec[0],
1068                                 fshift[0],
1069                                 &gid,
1070                                 x14[0],
1071                                 f14[0],
1072                                 chargeA,
1073                                 &eps,
1074                                 &krf,
1075                                 &crf,
1076                                 egcoul,
1077                                 typeA,
1078                                 &ntype,
1079                                 nbfp,
1080                                 egnb,
1081                                 &tabscale,
1082                                 tab,
1083                                 NULL,
1084                                 NULL,
1085                                 NULL,
1086                                 NULL,
1087                                 &nthreads,
1088                                 &count,
1089                                 (void *) &mtx,
1090                                 &outeriter,
1091                                 &inneriter,
1092                                 fr->adress_ex_forcecap,
1093                                 wf14);
1094                     } else {
1095                         nb_kernel330_adress_ex(&i1,
1096                                 &i0,
1097                                 j_index,
1098                                 &i1,
1099                                 &shift_f,
1100                                 fr->shift_vec[0],
1101                                 fshift[0],
1102                                 &gid,
1103                                 x14[0],
1104                                 f14[0],
1105                                 chargeA,
1106                                 &eps,
1107                                 &krf,
1108                                 &crf,
1109                                 egcoul,
1110                                 typeA,
1111                                 &ntype,
1112                                 nbfp,
1113                                 egnb,
1114                                 &tabscale,
1115                                 tab,
1116                                 NULL,
1117                                 NULL,
1118                                 NULL,
1119                                 NULL,
1120                                 &nthreads,
1121                                 &count,
1122                                 (void *) &mtx,
1123                                 &outeriter,
1124                                 &inneriter,
1125                                 fr->adress_ex_forcecap,
1126                                 wf14);
1127                     }
1128
1129                 }
1130             }
1131
1132         /* Add the forces */
1133         rvec_inc(f[ai],f14[0]);
1134         rvec_dec(f[aj],f14[0]);
1135
1136         if (g)
1137         {
1138             /* Correct the shift forces using the graph */
1139             ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1140             shift_vir = IVEC2IS(dt);
1141             rvec_inc(fshift[shift_vir],f14[0]);
1142             rvec_dec(fshift[CENTRAL],f14[0]);
1143         }
1144
1145             /* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
1146         }
1147     }
1148     return 0.0;
1149 }
1150
1151