src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel103_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_kernel103_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel103_adress_cg
  33  * Coulomb interaction:     Normal Coulomb
  34  * VdW interaction:         Not calculated
  35  * water optimization:      TIP4P - other atoms
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel103_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         rinvsq;
  79     real         jq;
  80     real         qq,vcoul,vctot;
  81     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  82     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  83     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
  84     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
  85     real         dx21,dy21,dz21,rsq21,rinv21;
  86     real         dx31,dy31,dz31,rsq31,rinv31;
  87     real         dx41,dy41,dz41,rsq41,rinv41;
  88     real         qH,qM;
  89     real         weight_cg1, weight_cg2, weight_product;
  90     real         hybscal;
  91
  92     nri              = *p_nri;
  93     ntype            = *p_ntype;
  94     nthreads         = *p_nthreads;
  95     facel            = *p_facel;
  96     krf              = *p_krf;
  97     crf              = *p_crf;
  98     tabscale         = *p_tabscale;
  99     ii               = iinr[0];
 100     qH               = facel*charge[ii+1];
 101     qM               = facel*charge[ii+3];
 102
 103     nouter           = 0;
 104     ninner           = 0;
 105
 106     do
 107     {
 108         #ifdef GMX_THREAD_SHM_FDECOMP
 109         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 110         nn0              = *count;
 111         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 112         *count           = nn1;
 113         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 114         if(nn1>nri) nn1=nri;
 115         #else
 116         nn0 = 0;
 117         nn1 = nri;
 118         #endif
 119
 120         for(n=nn0; (n<nn1); n++)
 121         {
 122             is3              = 3*shift[n];
 123             shX              = shiftvec[is3];
 124             shY              = shiftvec[is3+1];
 125             shZ              = shiftvec[is3+2];
 126             nj0              = jindex[n];
 127             nj1              = jindex[n+1];
 128             ii               = iinr[n];
 129             ii3              = 3*ii;
 130             ix2              = shX + pos[ii3+3];
 131             iy2              = shY + pos[ii3+4];
 132             iz2              = shZ + pos[ii3+5];
 133             ix3              = shX + pos[ii3+6];
 134             iy3              = shY + pos[ii3+7];
 135             iz3              = shZ + pos[ii3+8];
 136             ix4              = shX + pos[ii3+9];
 137             iy4              = shY + pos[ii3+10];
 138             iz4              = shZ + pos[ii3+11];
 139             weight_cg1       = wf[ii];
 140             vctot            = 0;
 141             fix2             = 0;
 142             fiy2             = 0;
 143             fiz2             = 0;
 144             fix3             = 0;
 145             fiy3             = 0;
 146             fiz3             = 0;
 147             fix4             = 0;
 148             fiy4             = 0;
 149             fiz4             = 0;
 150
 151             for(k=nj0; (k<nj1); k++)
 152             {
 153                 jnr              = jjnr[k];
 154                 weight_cg2       = wf[jnr];
 155                 weight_product   = weight_cg1*weight_cg2;
 156                 if (weight_product < ALMOST_ZERO) {
 157                        hybscal = 1.0;
 158                 }
 159                 else if (weight_product >= ALMOST_ONE)
 160                 {
 161                   /* force is zero, skip this molecule */
 162                        continue;
 163                 }
 164                 else
 165                 {
 166                    hybscal = 1.0 - weight_product;
 167                 }
 168                 j3               = 3*jnr;
 169                 jx1              = pos[j3+0];
 170                 jy1              = pos[j3+1];
 171                 jz1              = pos[j3+2];
 172                 dx21             = ix2 - jx1;
 173                 dy21             = iy2 - jy1;
 174                 dz21             = iz2 - jz1;
 175                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 176                 dx31             = ix3 - jx1;
 177                 dy31             = iy3 - jy1;
 178                 dz31             = iz3 - jz1;
 179                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 180                 dx41             = ix4 - jx1;
 181                 dy41             = iy4 - jy1;
 182                 dz41             = iz4 - jz1;
 183                 rsq41            = dx41*dx41+dy41*dy41+dz41*dz41;
 184                 rinv21           = 1.0/sqrt(rsq21);
 185                 rinv31           = 1.0/sqrt(rsq31);
 186                 rinv41           = 1.0/sqrt(rsq41);
 187                 jq               = charge[jnr+0];
 188                 qq               = qH*jq;
 189                 rinvsq           = rinv21*rinv21;
 190                 vcoul            = qq*rinv21;
 191                 vctot            = vctot+vcoul;
 192                 fscal            = (vcoul)*rinvsq;
 193                 fscal *= hybscal;
 194                 tx               = fscal*dx21;
 195                 ty               = fscal*dy21;
 196                 tz               = fscal*dz21;
 197                 fix2             = fix2 + tx;
 198                 fiy2             = fiy2 + ty;
 199                 fiz2             = fiz2 + tz;
 200                 fjx1             = faction[j3+0] - tx;
 201                 fjy1             = faction[j3+1] - ty;
 202                 fjz1             = faction[j3+2] - tz;
 203                 rinvsq           = rinv31*rinv31;
 204                 vcoul            = qq*rinv31;
 205                 vctot            = vctot+vcoul;
 206                 fscal            = (vcoul)*rinvsq;
 207                 fscal *= hybscal;
 208                 tx               = fscal*dx31;
 209                 ty               = fscal*dy31;
 210                 tz               = fscal*dz31;
 211                 fix3             = fix3 + tx;
 212                 fiy3             = fiy3 + ty;
 213                 fiz3             = fiz3 + tz;
 214                 fjx1             = fjx1 - tx;
 215                 fjy1             = fjy1 - ty;
 216                 fjz1             = fjz1 - tz;
 217                 qq               = qM*jq;
 218                 rinvsq           = rinv41*rinv41;
 219                 vcoul            = qq*rinv41;
 220                 vctot            = vctot+vcoul;
 221                 fscal            = (vcoul)*rinvsq;
 222                 fscal *= hybscal;
 223                 tx               = fscal*dx41;
 224                 ty               = fscal*dy41;
 225                 tz               = fscal*dz41;
 226                 fix4             = fix4 + tx;
 227                 fiy4             = fiy4 + ty;
 228                 fiz4             = fiz4 + tz;
 229                 faction[j3+0]    = fjx1 - tx;
 230                 faction[j3+1]    = fjy1 - ty;
 231                 faction[j3+2]    = fjz1 - tz;
 232             }
 233
 234             faction[ii3+3]   = faction[ii3+3] + fix2;
 235             faction[ii3+4]   = faction[ii3+4] + fiy2;
 236             faction[ii3+5]   = faction[ii3+5] + fiz2;
 237             faction[ii3+6]   = faction[ii3+6] + fix3;
 238             faction[ii3+7]   = faction[ii3+7] + fiy3;
 239             faction[ii3+8]   = faction[ii3+8] + fiz3;
 240             faction[ii3+9]   = faction[ii3+9] + fix4;
 241             faction[ii3+10]  = faction[ii3+10] + fiy4;
 242             faction[ii3+11]  = faction[ii3+11] + fiz4;
 243             fshift[is3]      = fshift[is3]+fix2+fix3+fix4;
 244             fshift[is3+1]    = fshift[is3+1]+fiy2+fiy3+fiy4;
 245             fshift[is3+2]    = fshift[is3+2]+fiz2+fiz3+fiz4;
 246             ggid             = gid[n];
 247             Vc[ggid]         = Vc[ggid] + vctot;
 248             ninner           = ninner + nj1 - nj0;
 249         }
 250
 251         nouter           = nouter + nn1 - nn0;
 252     }
 253     while (nn1<nri);
 254
 255     *outeriter       = nouter;
 256     *inneriter       = ninner;
 257 }
 258
 259
 260
 261
 262
 263 /*
 264  * Gromacs nonbonded kernel nb_kernel103_adress_ex
 265  * Coulomb interaction:     Normal Coulomb
 266  * VdW interaction:         Not calculated
 267  * water optimization:      TIP4P - other atoms
 268  * Calculate forces:        yes
 269  */
 270 void nb_kernel103_adress_ex(
 271                     int *           p_nri,
 272                     int *           iinr,
 273                     int *           jindex,
 274                     int *           jjnr,
 275                     int *           shift,
 276                     real *         shiftvec,
 277                     real *         fshift,
 278                     int *           gid,
 279                     real *         pos,
 280                     real *         faction,
 281                     real *         charge,
 282                     real *         p_facel,
 283                     real *         p_krf,
 284                     real *         p_crf,
 285                     real *         Vc,
 286                     int *           type,
 287                     int *           p_ntype,
 288                     real *         vdwparam,
 289                     real *         Vvdw,
 290                     real *         p_tabscale,
 291                     real *         VFtab,
 292                     real *         invsqrta,
 293                     real *         dvda,
 294                     real *         p_gbtabscale,
 295                     real *         GBtab,
 296                     int *           p_nthreads,
 297                     int *           count,
 298                     void *          mtx,
 299                     int *           outeriter,
 300                     int *           inneriter,
 301                     real           force_cap,
 302                     real *         wf)
 303 {
 304     int           nri,ntype,nthreads;
 305     real         facel,krf,crf,tabscale,gbtabscale;
 306     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 307     int           nn0,nn1,nouter,ninner;
 308     real         shX,shY,shZ;
 309     real         fscal,tx,ty,tz;
 310     real         rinvsq;
 311     real         jq;
 312     real         qq,vcoul,vctot;
 313     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 314     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 315     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
 316     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
 317     real         dx21,dy21,dz21,rsq21,rinv21;
 318     real         dx31,dy31,dz31,rsq31,rinv31;
 319     real         dx41,dy41,dz41,rsq41,rinv41;
 320     real         qH,qM;
 321     real         weight_cg1, weight_cg2, weight_product;
 322     real         hybscal;
 323
 324     nri              = *p_nri;
 325     ntype            = *p_ntype;
 326     nthreads         = *p_nthreads;
 327     facel            = *p_facel;
 328     krf              = *p_krf;
 329     crf              = *p_crf;
 330     tabscale         = *p_tabscale;
 331     ii               = iinr[0];
 332     qH               = facel*charge[ii+1];
 333     qM               = facel*charge[ii+3];
 334
 335     nouter           = 0;
 336     ninner           = 0;
 337
 338     do
 339     {
 340         #ifdef GMX_THREAD_SHM_FDECOMP
 341         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 342         nn0              = *count;
 343         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 344         *count           = nn1;
 345         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 346         if(nn1>nri) nn1=nri;
 347         #else
 348         nn0 = 0;
 349         nn1 = nri;
 350         #endif
 351
 352         for(n=nn0; (n<nn1); n++)
 353         {
 354             is3              = 3*shift[n];
 355             shX              = shiftvec[is3];
 356             shY              = shiftvec[is3+1];
 357             shZ              = shiftvec[is3+2];
 358             nj0              = jindex[n];
 359             nj1              = jindex[n+1];
 360             ii               = iinr[n];
 361             ii3              = 3*ii;
 362             ix2              = shX + pos[ii3+3];
 363             iy2              = shY + pos[ii3+4];
 364             iz2              = shZ + pos[ii3+5];
 365             ix3              = shX + pos[ii3+6];
 366             iy3              = shY + pos[ii3+7];
 367             iz3              = shZ + pos[ii3+8];
 368             ix4              = shX + pos[ii3+9];
 369             iy4              = shY + pos[ii3+10];
 370             iz4              = shZ + pos[ii3+11];
 371             weight_cg1       = wf[ii];
 372             vctot            = 0;
 373             fix2             = 0;
 374             fiy2             = 0;
 375             fiz2             = 0;
 376             fix3             = 0;
 377             fiy3             = 0;
 378             fiz3             = 0;
 379             fix4             = 0;
 380             fiy4             = 0;
 381             fiz4             = 0;
 382
 383             for(k=nj0; (k<nj1); k++)
 384             {
 385                 jnr              = jjnr[k];
 386                 weight_cg2       = wf[jnr];
 387                 weight_product   = weight_cg1*weight_cg2;
 388                 if (weight_product < ALMOST_ZERO) {
 389                 /* force is zero, skip this molecule */
 390                  continue;
 391                 }
 392                 else if (weight_product >= ALMOST_ONE)
 393                 {
 394                        hybscal = 1.0;
 395                 }
 396                 else
 397                 {
 398                    hybscal = weight_product;
 399                 }
 400                 j3               = 3*jnr;
 401                 jx1              = pos[j3+0];
 402                 jy1              = pos[j3+1];
 403                 jz1              = pos[j3+2];
 404                 dx21             = ix2 - jx1;
 405                 dy21             = iy2 - jy1;
 406                 dz21             = iz2 - jz1;
 407                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 408                 dx31             = ix3 - jx1;
 409                 dy31             = iy3 - jy1;
 410                 dz31             = iz3 - jz1;
 411                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 412                 dx41             = ix4 - jx1;
 413                 dy41             = iy4 - jy1;
 414                 dz41             = iz4 - jz1;
 415                 rsq41            = dx41*dx41+dy41*dy41+dz41*dz41;
 416                 rinv21           = 1.0/sqrt(rsq21);
 417                 rinv31           = 1.0/sqrt(rsq31);
 418                 rinv41           = 1.0/sqrt(rsq41);
 419                 jq               = charge[jnr+0];
 420                 qq               = qH*jq;
 421                 rinvsq           = rinv21*rinv21;
 422                 vcoul            = qq*rinv21;
 423                 vctot            = vctot+vcoul;
 424                 fscal            = (vcoul)*rinvsq;
 425                 fscal *= hybscal;
 426                 tx               = fscal*dx21;
 427                 ty               = fscal*dy21;
 428                 tz               = fscal*dz21;
 429                 fix2             = fix2 + tx;
 430                 fiy2             = fiy2 + ty;
 431                 fiz2             = fiz2 + tz;
 432                 fjx1             = faction[j3+0] - tx;
 433                 fjy1             = faction[j3+1] - ty;
 434                 fjz1             = faction[j3+2] - tz;
 435                 rinvsq           = rinv31*rinv31;
 436                 vcoul            = qq*rinv31;
 437                 vctot            = vctot+vcoul;
 438                 fscal            = (vcoul)*rinvsq;
 439                 fscal *= hybscal;
 440                 tx               = fscal*dx31;
 441                 ty               = fscal*dy31;
 442                 tz               = fscal*dz31;
 443                 fix3             = fix3 + tx;
 444                 fiy3             = fiy3 + ty;
 445                 fiz3             = fiz3 + tz;
 446                 fjx1             = fjx1 - tx;
 447                 fjy1             = fjy1 - ty;
 448                 fjz1             = fjz1 - tz;
 449                 qq               = qM*jq;
 450                 rinvsq           = rinv41*rinv41;
 451                 vcoul            = qq*rinv41;
 452                 vctot            = vctot+vcoul;
 453                 fscal            = (vcoul)*rinvsq;
 454                 fscal *= hybscal;
 455                 tx               = fscal*dx41;
 456                 ty               = fscal*dy41;
 457                 tz               = fscal*dz41;
 458                 fix4             = fix4 + tx;
 459                 fiy4             = fiy4 + ty;
 460                 fiz4             = fiz4 + tz;
 461                 faction[j3+0]    = fjx1 - tx;
 462                 faction[j3+1]    = fjy1 - ty;
 463                 faction[j3+2]    = fjz1 - tz;
 464             }
 465
 466             faction[ii3+3]   = faction[ii3+3] + fix2;
 467             faction[ii3+4]   = faction[ii3+4] + fiy2;
 468             faction[ii3+5]   = faction[ii3+5] + fiz2;
 469             faction[ii3+6]   = faction[ii3+6] + fix3;
 470             faction[ii3+7]   = faction[ii3+7] + fiy3;
 471             faction[ii3+8]   = faction[ii3+8] + fiz3;
 472             faction[ii3+9]   = faction[ii3+9] + fix4;
 473             faction[ii3+10]  = faction[ii3+10] + fiy4;
 474             faction[ii3+11]  = faction[ii3+11] + fiz4;
 475             fshift[is3]      = fshift[is3]+fix2+fix3+fix4;
 476             fshift[is3+1]    = fshift[is3+1]+fiy2+fiy3+fiy4;
 477             fshift[is3+2]    = fshift[is3+2]+fiz2+fiz3+fiz4;
 478             ggid             = gid[n];
 479             Vc[ggid]         = Vc[ggid] + vctot;
 480             ninner           = ninner + nj1 - nj0;
 481         }
 482
 483         nouter           = nouter + nn1 - nn0;
 484     }
 485     while (nn1<nri);
 486
 487     *outeriter       = nouter;
 488     *inneriter       = ninner;
 489 }
 490
 491