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