src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel301_c_adress.c

   1 /*
   2  * Copyright (c) Erik Lindahl, David van der Spoel 2003
   3  *
   4  * This file is generated automatically at compile time
   5  * by the program mknb in the Gromacs distribution.
   6  *
   7  * Options used when generation this file:
   8  * Language:         c
   9  * Precision:        single
  10  * Threads:          yes
  11  * Software invsqrt: no
  12  * PowerPC invsqrt:  no
  13  * Prefetch forces:  no
  14  * Adress kernel:  yes
  15  * Comments:         no
  16  */
  17 #ifdef HAVE_CONFIG_H
  18 #include<config.h>
  19 #endif
  20 #ifdef GMX_THREAD_SHM_FDECOMP
  21 #include<thread_mpi.h>
  22 #endif
  23 #define ALMOST_ZERO 1e-30
  24 #define ALMOST_ONE 1-(1e-30)
  25 #include<math.h>
  26
  27 #include "nb_kernel301_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel301_adress_cg
  33  * Coulomb interaction:     Tabulated
  34  * VdW interaction:         Not calculated
  35  * water optimization:      SPC/TIP3P - other atoms
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel301_adress_cg(
  39                     int *           p_nri,
  40                     int *           iinr,
  41                     int *           jindex,
  42                     int *           jjnr,
  43                     int *           shift,
  44                     real *         shiftvec,
  45                     real *         fshift,
  46                     int *           gid,
  47                     real *         pos,
  48                     real *         faction,
  49                     real *         charge,
  50                     real *         p_facel,
  51                     real *         p_krf,
  52                     real *         p_crf,
  53                     real *         Vc,
  54                     int *           type,
  55                     int *           p_ntype,
  56                     real *         vdwparam,
  57                     real *         Vvdw,
  58                     real *         p_tabscale,
  59                     real *         VFtab,
  60                     real *         invsqrta,
  61                     real *         dvda,
  62                     real *         p_gbtabscale,
  63                     real *         GBtab,
  64                     int *           p_nthreads,
  65                     int *           count,
  66                     void *          mtx,
  67                     int *           outeriter,
  68                     int *           inneriter,
  69                     real           force_cap,
  70                     real *         wf)
  71 {
  72     int           nri,ntype,nthreads;
  73     real         facel,krf,crf,tabscale,gbtabscale;
  74     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
  75     int           nn0,nn1,nouter,ninner;
  76     real         shX,shY,shZ;
  77     real         fscal,tx,ty,tz;
  78     real         jq;
  79     real         qq,vcoul,vctot;
  80     real         r,rt,eps,eps2;
  81     int           n0,nnn;
  82     real         Y,F,Geps,Heps2,Fp,VV;
  83     real         FF;
  84     real         fijC;
  85     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  86     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  87     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  88     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
  89     real         dx11,dy11,dz11,rsq11,rinv11;
  90     real         dx21,dy21,dz21,rsq21,rinv21;
  91     real         dx31,dy31,dz31,rsq31,rinv31;
  92     real         qO,qH;
  93     real         weight_cg1, weight_cg2, weight_product;
  94     real         hybscal;
  95
  96     nri              = *p_nri;
  97     ntype            = *p_ntype;
  98     nthreads         = *p_nthreads;
  99     facel            = *p_facel;
 100     krf              = *p_krf;
 101     crf              = *p_crf;
 102     tabscale         = *p_tabscale;
 103     ii               = iinr[0];
 104     qO               = facel*charge[ii];
 105     qH               = facel*charge[ii+1];
 106
 107     nouter           = 0;
 108     ninner           = 0;
 109
 110     do
 111     {
 112         #ifdef GMX_THREAD_SHM_FDECOMP
 113         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 114         nn0              = *count;
 115         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 116         *count           = nn1;
 117         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 118         if(nn1>nri) nn1=nri;
 119         #else
 120         nn0 = 0;
 121         nn1 = nri;
 122         #endif
 123
 124         for(n=nn0; (n<nn1); n++)
 125         {
 126             is3              = 3*shift[n];
 127             shX              = shiftvec[is3];
 128             shY              = shiftvec[is3+1];
 129             shZ              = shiftvec[is3+2];
 130             nj0              = jindex[n];
 131             nj1              = jindex[n+1];
 132             ii               = iinr[n];
 133             ii3              = 3*ii;
 134             ix1              = shX + pos[ii3+0];
 135             iy1              = shY + pos[ii3+1];
 136             iz1              = shZ + pos[ii3+2];
 137             ix2              = shX + pos[ii3+3];
 138             iy2              = shY + pos[ii3+4];
 139             iz2              = shZ + pos[ii3+5];
 140             ix3              = shX + pos[ii3+6];
 141             iy3              = shY + pos[ii3+7];
 142             iz3              = shZ + pos[ii3+8];
 143             weight_cg1       = wf[ii];
 144             vctot            = 0;
 145             fix1             = 0;
 146             fiy1             = 0;
 147             fiz1             = 0;
 148             fix2             = 0;
 149             fiy2             = 0;
 150             fiz2             = 0;
 151             fix3             = 0;
 152             fiy3             = 0;
 153             fiz3             = 0;
 154
 155             for(k=nj0; (k<nj1); k++)
 156             {
 157                 jnr              = jjnr[k];
 158                 weight_cg2       = wf[jnr];
 159                 weight_product   = weight_cg1*weight_cg2;
 160                 if (weight_product < ALMOST_ZERO) {
 161                        hybscal = 1.0;
 162                 }
 163                 else if (weight_product >= ALMOST_ONE)
 164                 {
 165                   /* force is zero, skip this molecule */
 166                        continue;
 167                 }
 168                 else
 169                 {
 170                    hybscal = 1.0 - weight_product;
 171                 }
 172                 j3               = 3*jnr;
 173                 jx1              = pos[j3+0];
 174                 jy1              = pos[j3+1];
 175                 jz1              = pos[j3+2];
 176                 dx11             = ix1 - jx1;
 177                 dy11             = iy1 - jy1;
 178                 dz11             = iz1 - jz1;
 179                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 180                 dx21             = ix2 - jx1;
 181                 dy21             = iy2 - jy1;
 182                 dz21             = iz2 - jz1;
 183                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 184                 dx31             = ix3 - jx1;
 185                 dy31             = iy3 - jy1;
 186                 dz31             = iz3 - jz1;
 187                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 188                 rinv11           = 1.0/sqrt(rsq11);
 189                 rinv21           = 1.0/sqrt(rsq21);
 190                 rinv31           = 1.0/sqrt(rsq31);
 191                 jq               = charge[jnr+0];
 192                 qq               = qO*jq;
 193                 r                = rsq11*rinv11;
 194                 rt               = r*tabscale;
 195                 n0               = rt;
 196                 eps              = rt-n0;
 197                 eps2             = eps*eps;
 198                 nnn              = 4*n0;
 199                 Y                = VFtab[nnn];
 200                 F                = VFtab[nnn+1];
 201                 Geps             = eps*VFtab[nnn+2];
 202                 Heps2            = eps2*VFtab[nnn+3];
 203                 Fp               = F+Geps+Heps2;
 204                 VV               = Y+eps*Fp;
 205                 FF               = Fp+Geps+2.0*Heps2;
 206                 vcoul            = qq*VV;
 207                 fijC             = qq*FF;
 208                 vctot            = vctot + vcoul;
 209                 fscal            = -((fijC)*tabscale)*rinv11;
 210                 fscal *= hybscal;
 211                 tx               = fscal*dx11;
 212                 ty               = fscal*dy11;
 213                 tz               = fscal*dz11;
 214                 fix1             = fix1 + tx;
 215                 fiy1             = fiy1 + ty;
 216                 fiz1             = fiz1 + tz;
 217                 fjx1             = faction[j3+0] - tx;
 218                 fjy1             = faction[j3+1] - ty;
 219                 fjz1             = faction[j3+2] - tz;
 220                 qq               = qH*jq;
 221                 r                = rsq21*rinv21;
 222                 rt               = r*tabscale;
 223                 n0               = rt;
 224                 eps              = rt-n0;
 225                 eps2             = eps*eps;
 226                 nnn              = 4*n0;
 227                 Y                = VFtab[nnn];
 228                 F                = VFtab[nnn+1];
 229                 Geps             = eps*VFtab[nnn+2];
 230                 Heps2            = eps2*VFtab[nnn+3];
 231                 Fp               = F+Geps+Heps2;
 232                 VV               = Y+eps*Fp;
 233                 FF               = Fp+Geps+2.0*Heps2;
 234                 vcoul            = qq*VV;
 235                 fijC             = qq*FF;
 236                 vctot            = vctot + vcoul;
 237                 fscal            = -((fijC)*tabscale)*rinv21;
 238                 fscal *= hybscal;
 239                 tx               = fscal*dx21;
 240                 ty               = fscal*dy21;
 241                 tz               = fscal*dz21;
 242                 fix2             = fix2 + tx;
 243                 fiy2             = fiy2 + ty;
 244                 fiz2             = fiz2 + tz;
 245                 fjx1             = fjx1 - tx;
 246                 fjy1             = fjy1 - ty;
 247                 fjz1             = fjz1 - tz;
 248                 r                = rsq31*rinv31;
 249                 rt               = r*tabscale;
 250                 n0               = rt;
 251                 eps              = rt-n0;
 252                 eps2             = eps*eps;
 253                 nnn              = 4*n0;
 254                 Y                = VFtab[nnn];
 255                 F                = VFtab[nnn+1];
 256                 Geps             = eps*VFtab[nnn+2];
 257                 Heps2            = eps2*VFtab[nnn+3];
 258                 Fp               = F+Geps+Heps2;
 259                 VV               = Y+eps*Fp;
 260                 FF               = Fp+Geps+2.0*Heps2;
 261                 vcoul            = qq*VV;
 262                 fijC             = qq*FF;
 263                 vctot            = vctot + vcoul;
 264                 fscal            = -((fijC)*tabscale)*rinv31;
 265                 fscal *= hybscal;
 266                 tx               = fscal*dx31;
 267                 ty               = fscal*dy31;
 268                 tz               = fscal*dz31;
 269                 fix3             = fix3 + tx;
 270                 fiy3             = fiy3 + ty;
 271                 fiz3             = fiz3 + tz;
 272                 faction[j3+0]    = fjx1 - tx;
 273                 faction[j3+1]    = fjy1 - ty;
 274                 faction[j3+2]    = fjz1 - tz;
 275             }
 276
 277             faction[ii3+0]   = faction[ii3+0] + fix1;
 278             faction[ii3+1]   = faction[ii3+1] + fiy1;
 279             faction[ii3+2]   = faction[ii3+2] + fiz1;
 280             faction[ii3+3]   = faction[ii3+3] + fix2;
 281             faction[ii3+4]   = faction[ii3+4] + fiy2;
 282             faction[ii3+5]   = faction[ii3+5] + fiz2;
 283             faction[ii3+6]   = faction[ii3+6] + fix3;
 284             faction[ii3+7]   = faction[ii3+7] + fiy3;
 285             faction[ii3+8]   = faction[ii3+8] + fiz3;
 286             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 287             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 288             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 289             ggid             = gid[n];
 290             Vc[ggid]         = Vc[ggid] + vctot;
 291             ninner           = ninner + nj1 - nj0;
 292         }
 293
 294         nouter           = nouter + nn1 - nn0;
 295     }
 296     while (nn1<nri);
 297
 298     *outeriter       = nouter;
 299     *inneriter       = ninner;
 300 }
 301
 302
 303
 304
 305
 306 /*
 307  * Gromacs nonbonded kernel nb_kernel301_adress_ex
 308  * Coulomb interaction:     Tabulated
 309  * VdW interaction:         Not calculated
 310  * water optimization:      SPC/TIP3P - other atoms
 311  * Calculate forces:        yes
 312  */
 313 void nb_kernel301_adress_ex(
 314                     int *           p_nri,
 315                     int *           iinr,
 316                     int *           jindex,
 317                     int *           jjnr,
 318                     int *           shift,
 319                     real *         shiftvec,
 320                     real *         fshift,
 321                     int *           gid,
 322                     real *         pos,
 323                     real *         faction,
 324                     real *         charge,
 325                     real *         p_facel,
 326                     real *         p_krf,
 327                     real *         p_crf,
 328                     real *         Vc,
 329                     int *           type,
 330                     int *           p_ntype,
 331                     real *         vdwparam,
 332                     real *         Vvdw,
 333                     real *         p_tabscale,
 334                     real *         VFtab,
 335                     real *         invsqrta,
 336                     real *         dvda,
 337                     real *         p_gbtabscale,
 338                     real *         GBtab,
 339                     int *           p_nthreads,
 340                     int *           count,
 341                     void *          mtx,
 342                     int *           outeriter,
 343                     int *           inneriter,
 344                     real           force_cap,
 345                     real *         wf)
 346 {
 347     int           nri,ntype,nthreads;
 348     real         facel,krf,crf,tabscale,gbtabscale;
 349     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 350     int           nn0,nn1,nouter,ninner;
 351     real         shX,shY,shZ;
 352     real         fscal,tx,ty,tz;
 353     real         jq;
 354     real         qq,vcoul,vctot;
 355     real         r,rt,eps,eps2;
 356     int           n0,nnn;
 357     real         Y,F,Geps,Heps2,Fp,VV;
 358     real         FF;
 359     real         fijC;
 360     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 361     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 362     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 363     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
 364     real         dx11,dy11,dz11,rsq11,rinv11;
 365     real         dx21,dy21,dz21,rsq21,rinv21;
 366     real         dx31,dy31,dz31,rsq31,rinv31;
 367     real         qO,qH;
 368     real         weight_cg1, weight_cg2, weight_product;
 369     real         hybscal;
 370
 371     nri              = *p_nri;
 372     ntype            = *p_ntype;
 373     nthreads         = *p_nthreads;
 374     facel            = *p_facel;
 375     krf              = *p_krf;
 376     crf              = *p_crf;
 377     tabscale         = *p_tabscale;
 378     ii               = iinr[0];
 379     qO               = facel*charge[ii];
 380     qH               = facel*charge[ii+1];
 381
 382     nouter           = 0;
 383     ninner           = 0;
 384
 385     do
 386     {
 387         #ifdef GMX_THREAD_SHM_FDECOMP
 388         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 389         nn0              = *count;
 390         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 391         *count           = nn1;
 392         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 393         if(nn1>nri) nn1=nri;
 394         #else
 395         nn0 = 0;
 396         nn1 = nri;
 397         #endif
 398
 399         for(n=nn0; (n<nn1); n++)
 400         {
 401             is3              = 3*shift[n];
 402             shX              = shiftvec[is3];
 403             shY              = shiftvec[is3+1];
 404             shZ              = shiftvec[is3+2];
 405             nj0              = jindex[n];
 406             nj1              = jindex[n+1];
 407             ii               = iinr[n];
 408             ii3              = 3*ii;
 409             ix1              = shX + pos[ii3+0];
 410             iy1              = shY + pos[ii3+1];
 411             iz1              = shZ + pos[ii3+2];
 412             ix2              = shX + pos[ii3+3];
 413             iy2              = shY + pos[ii3+4];
 414             iz2              = shZ + pos[ii3+5];
 415             ix3              = shX + pos[ii3+6];
 416             iy3              = shY + pos[ii3+7];
 417             iz3              = shZ + pos[ii3+8];
 418             weight_cg1       = wf[ii];
 419             vctot            = 0;
 420             fix1             = 0;
 421             fiy1             = 0;
 422             fiz1             = 0;
 423             fix2             = 0;
 424             fiy2             = 0;
 425             fiz2             = 0;
 426             fix3             = 0;
 427             fiy3             = 0;
 428             fiz3             = 0;
 429
 430             for(k=nj0; (k<nj1); k++)
 431             {
 432                 jnr              = jjnr[k];
 433                 weight_cg2       = wf[jnr];
 434                 weight_product   = weight_cg1*weight_cg2;
 435                 if (weight_product < ALMOST_ZERO) {
 436                 /* force is zero, skip this molecule */
 437                  continue;
 438                 }
 439                 else if (weight_product >= ALMOST_ONE)
 440                 {
 441                        hybscal = 1.0;
 442                 }
 443                 else
 444                 {
 445                    hybscal = weight_product;
 446                 }
 447                 j3               = 3*jnr;
 448                 jx1              = pos[j3+0];
 449                 jy1              = pos[j3+1];
 450                 jz1              = pos[j3+2];
 451                 dx11             = ix1 - jx1;
 452                 dy11             = iy1 - jy1;
 453                 dz11             = iz1 - jz1;
 454                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 455                 dx21             = ix2 - jx1;
 456                 dy21             = iy2 - jy1;
 457                 dz21             = iz2 - jz1;
 458                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 459                 dx31             = ix3 - jx1;
 460                 dy31             = iy3 - jy1;
 461                 dz31             = iz3 - jz1;
 462                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 463                 rinv11           = 1.0/sqrt(rsq11);
 464                 rinv21           = 1.0/sqrt(rsq21);
 465                 rinv31           = 1.0/sqrt(rsq31);
 466                 jq               = charge[jnr+0];
 467                 qq               = qO*jq;
 468                 r                = rsq11*rinv11;
 469                 rt               = r*tabscale;
 470                 n0               = rt;
 471                 eps              = rt-n0;
 472                 eps2             = eps*eps;
 473                 nnn              = 4*n0;
 474                 Y                = VFtab[nnn];
 475                 F                = VFtab[nnn+1];
 476                 Geps             = eps*VFtab[nnn+2];
 477                 Heps2            = eps2*VFtab[nnn+3];
 478                 Fp               = F+Geps+Heps2;
 479                 VV               = Y+eps*Fp;
 480                 FF               = Fp+Geps+2.0*Heps2;
 481                 vcoul            = qq*VV;
 482                 fijC             = qq*FF;
 483                 vctot            = vctot + vcoul;
 484                 fscal            = -((fijC)*tabscale)*rinv11;
 485                 fscal *= hybscal;
 486                 tx               = fscal*dx11;
 487                 ty               = fscal*dy11;
 488                 tz               = fscal*dz11;
 489                 fix1             = fix1 + tx;
 490                 fiy1             = fiy1 + ty;
 491                 fiz1             = fiz1 + tz;
 492                 fjx1             = faction[j3+0] - tx;
 493                 fjy1             = faction[j3+1] - ty;
 494                 fjz1             = faction[j3+2] - tz;
 495                 qq               = qH*jq;
 496                 r                = rsq21*rinv21;
 497                 rt               = r*tabscale;
 498                 n0               = rt;
 499                 eps              = rt-n0;
 500                 eps2             = eps*eps;
 501                 nnn              = 4*n0;
 502                 Y                = VFtab[nnn];
 503                 F                = VFtab[nnn+1];
 504                 Geps             = eps*VFtab[nnn+2];
 505                 Heps2            = eps2*VFtab[nnn+3];
 506                 Fp               = F+Geps+Heps2;
 507                 VV               = Y+eps*Fp;
 508                 FF               = Fp+Geps+2.0*Heps2;
 509                 vcoul            = qq*VV;
 510                 fijC             = qq*FF;
 511                 vctot            = vctot + vcoul;
 512                 fscal            = -((fijC)*tabscale)*rinv21;
 513                 fscal *= hybscal;
 514                 tx               = fscal*dx21;
 515                 ty               = fscal*dy21;
 516                 tz               = fscal*dz21;
 517                 fix2             = fix2 + tx;
 518                 fiy2             = fiy2 + ty;
 519                 fiz2             = fiz2 + tz;
 520                 fjx1             = fjx1 - tx;
 521                 fjy1             = fjy1 - ty;
 522                 fjz1             = fjz1 - tz;
 523                 r                = rsq31*rinv31;
 524                 rt               = r*tabscale;
 525                 n0               = rt;
 526                 eps              = rt-n0;
 527                 eps2             = eps*eps;
 528                 nnn              = 4*n0;
 529                 Y                = VFtab[nnn];
 530                 F                = VFtab[nnn+1];
 531                 Geps             = eps*VFtab[nnn+2];
 532                 Heps2            = eps2*VFtab[nnn+3];
 533                 Fp               = F+Geps+Heps2;
 534                 VV               = Y+eps*Fp;
 535                 FF               = Fp+Geps+2.0*Heps2;
 536                 vcoul            = qq*VV;
 537                 fijC             = qq*FF;
 538                 vctot            = vctot + vcoul;
 539                 fscal            = -((fijC)*tabscale)*rinv31;
 540                 fscal *= hybscal;
 541                 tx               = fscal*dx31;
 542                 ty               = fscal*dy31;
 543                 tz               = fscal*dz31;
 544                 fix3             = fix3 + tx;
 545                 fiy3             = fiy3 + ty;
 546                 fiz3             = fiz3 + tz;
 547                 faction[j3+0]    = fjx1 - tx;
 548                 faction[j3+1]    = fjy1 - ty;
 549                 faction[j3+2]    = fjz1 - tz;
 550             }
 551
 552             faction[ii3+0]   = faction[ii3+0] + fix1;
 553             faction[ii3+1]   = faction[ii3+1] + fiy1;
 554             faction[ii3+2]   = faction[ii3+2] + fiz1;
 555             faction[ii3+3]   = faction[ii3+3] + fix2;
 556             faction[ii3+4]   = faction[ii3+4] + fiy2;
 557             faction[ii3+5]   = faction[ii3+5] + fiz2;
 558             faction[ii3+6]   = faction[ii3+6] + fix3;
 559             faction[ii3+7]   = faction[ii3+7] + fiy3;
 560             faction[ii3+8]   = faction[ii3+8] + fiz3;
 561             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 562             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 563             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 564             ggid             = gid[n];
 565             Vc[ggid]         = Vc[ggid] + vctot;
 566             ninner           = ninner + nj1 - nj0;
 567         }
 568
 569         nouter           = nouter + nn1 - nn0;
 570     }
 571     while (nn1<nri);
 572
 573     *outeriter       = nouter;
 574     *inneriter       = ninner;
 575 }
 576
 577