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