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