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