src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel400_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_kernel400_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel400_adress_cg
  33  * Coulomb interaction:     Generalized-Born
  34  * VdW interaction:         Not calculated
  35  * water optimization:      No
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel400_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         iq;
  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         isai,isaj,isaprod,gbscale,vgb;
  86     real         dvdasum,dvdatmp,dvdaj,fgb;
  87     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  88     real         jx1,jy1,jz1;
  89     real         dx11,dy11,dz11,rsq11,rinv11;
  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     gbtabscale       = *p_gbtabscale;
 101     nouter           = 0;
 102     ninner           = 0;
 103
 104     do
 105     {
 106         #ifdef GMX_THREAD_SHM_FDECOMP
 107         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 108         nn0              = *count;
 109         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 110         *count           = nn1;
 111         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 112         if(nn1>nri) nn1=nri;
 113         #else
 114         nn0 = 0;
 115         nn1 = nri;
 116         #endif
 117
 118         for(n=nn0; (n<nn1); n++)
 119         {
 120             is3              = 3*shift[n];
 121             shX              = shiftvec[is3];
 122             shY              = shiftvec[is3+1];
 123             shZ              = shiftvec[is3+2];
 124             nj0              = jindex[n];
 125             nj1              = jindex[n+1];
 126             ii               = iinr[n];
 127             ii3              = 3*ii;
 128             ix1              = shX + pos[ii3+0];
 129             iy1              = shY + pos[ii3+1];
 130             iz1              = shZ + pos[ii3+2];
 131             iq               = facel*charge[ii];
 132             isai             = invsqrta[ii];
 133             weight_cg1       = wf[ii];
 134             vctot            = 0;
 135             dvdasum          = 0;
 136             fix1             = 0;
 137             fiy1             = 0;
 138             fiz1             = 0;
 139
 140             for(k=nj0; (k<nj1); k++)
 141             {
 142                 jnr              = jjnr[k];
 143                 weight_cg2       = wf[jnr];
 144                 weight_product   = weight_cg1*weight_cg2;
 145                 if (weight_product < ALMOST_ZERO) {
 146                        hybscal = 1.0;
 147                 }
 148                 else if (weight_product >= ALMOST_ONE)
 149                 {
 150                   /* force is zero, skip this molecule */
 151                        continue;
 152                 }
 153                 else
 154                 {
 155                    hybscal = 1.0 - weight_product;
 156                 }
 157                 j3               = 3*jnr;
 158                 jx1              = pos[j3+0];
 159                 jy1              = pos[j3+1];
 160                 jz1              = pos[j3+2];
 161                 dx11             = ix1 - jx1;
 162                 dy11             = iy1 - jy1;
 163                 dz11             = iz1 - jz1;
 164                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 165                 rinv11           = 1.0/sqrt(rsq11);
 166                 isaj             = invsqrta[jnr];
 167                 isaprod          = isai*isaj;
 168                 qq               = iq*charge[jnr];
 169                 vcoul            = qq*rinv11;
 170                 fscal            = vcoul*rinv11;
 171                 qq               = isaprod*(-qq);
 172                 gbscale          = isaprod*gbtabscale;
 173                 dvdaj            = dvda[jnr];
 174                 r                = rsq11*rinv11;
 175                 rt               = r*gbscale;
 176                 n0               = rt;
 177                 eps              = rt-n0;
 178                 eps2             = eps*eps;
 179                 nnn              = 4*n0;
 180                 Y                = GBtab[nnn];
 181                 F                = GBtab[nnn+1];
 182                 Geps             = eps*GBtab[nnn+2];
 183                 Heps2            = eps2*GBtab[nnn+3];
 184                 Fp               = F+Geps+Heps2;
 185                 VV               = Y+eps*Fp;
 186                 FF               = Fp+Geps+2.0*Heps2;
 187                 vgb              = qq*VV;
 188                 fijC             = qq*FF*gbscale;
 189                 dvdatmp          = -0.5*(vgb+fijC*r);
 190                 dvdasum          = dvdasum + dvdatmp;
 191                 dvda[jnr]        = dvdaj+dvdatmp*isaj*isaj;
 192                 vctot            = vctot + vcoul;
 193                 fscal            = -(fijC-fscal)*rinv11;
 194                 fscal *= hybscal;
 195                 tx               = fscal*dx11;
 196                 ty               = fscal*dy11;
 197                 tz               = fscal*dz11;
 198                 fix1             = fix1 + tx;
 199                 fiy1             = fiy1 + ty;
 200                 fiz1             = fiz1 + tz;
 201                 faction[j3+0]    = faction[j3+0] - tx;
 202                 faction[j3+1]    = faction[j3+1] - ty;
 203                 faction[j3+2]    = faction[j3+2] - tz;
 204             }
 205
 206             faction[ii3+0]   = faction[ii3+0] + fix1;
 207             faction[ii3+1]   = faction[ii3+1] + fiy1;
 208             faction[ii3+2]   = faction[ii3+2] + fiz1;
 209             fshift[is3]      = fshift[is3]+fix1;
 210             fshift[is3+1]    = fshift[is3+1]+fiy1;
 211             fshift[is3+2]    = fshift[is3+2]+fiz1;
 212             ggid             = gid[n];
 213             Vc[ggid]         = Vc[ggid] + vctot;
 214             dvda[ii]         = dvda[ii] + dvdasum*isai*isai;
 215             ninner           = ninner + nj1 - nj0;
 216         }
 217
 218         nouter           = nouter + nn1 - nn0;
 219     }
 220     while (nn1<nri);
 221
 222     *outeriter       = nouter;
 223     *inneriter       = ninner;
 224 }
 225
 226
 227
 228
 229
 230 /*
 231  * Gromacs nonbonded kernel nb_kernel400_adress_ex
 232  * Coulomb interaction:     Generalized-Born
 233  * VdW interaction:         Not calculated
 234  * water optimization:      No
 235  * Calculate forces:        yes
 236  */
 237 void nb_kernel400_adress_ex(
 238                     int *           p_nri,
 239                     int *           iinr,
 240                     int *           jindex,
 241                     int *           jjnr,
 242                     int *           shift,
 243                     real *         shiftvec,
 244                     real *         fshift,
 245                     int *           gid,
 246                     real *         pos,
 247                     real *         faction,
 248                     real *         charge,
 249                     real *         p_facel,
 250                     real *         p_krf,
 251                     real *         p_crf,
 252                     real *         Vc,
 253                     int *           type,
 254                     int *           p_ntype,
 255                     real *         vdwparam,
 256                     real *         Vvdw,
 257                     real *         p_tabscale,
 258                     real *         VFtab,
 259                     real *         invsqrta,
 260                     real *         dvda,
 261                     real *         p_gbtabscale,
 262                     real *         GBtab,
 263                     int *           p_nthreads,
 264                     int *           count,
 265                     void *          mtx,
 266                     int *           outeriter,
 267                     int *           inneriter,
 268                     real           force_cap,
 269                     real *         wf)
 270 {
 271     int           nri,ntype,nthreads;
 272     real         facel,krf,crf,tabscale,gbtabscale;
 273     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 274     int           nn0,nn1,nouter,ninner;
 275     real         shX,shY,shZ;
 276     real         fscal,tx,ty,tz;
 277     real         iq;
 278     real         qq,vcoul,vctot;
 279     real         r,rt,eps,eps2;
 280     int           n0,nnn;
 281     real         Y,F,Geps,Heps2,Fp,VV;
 282     real         FF;
 283     real         fijC;
 284     real         isai,isaj,isaprod,gbscale,vgb;
 285     real         dvdasum,dvdatmp,dvdaj,fgb;
 286     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 287     real         jx1,jy1,jz1;
 288     real         dx11,dy11,dz11,rsq11,rinv11;
 289     real         weight_cg1, weight_cg2, weight_product;
 290     real         hybscal;
 291
 292     nri              = *p_nri;
 293     ntype            = *p_ntype;
 294     nthreads         = *p_nthreads;
 295     facel            = *p_facel;
 296     krf              = *p_krf;
 297     crf              = *p_crf;
 298     tabscale         = *p_tabscale;
 299     gbtabscale       = *p_gbtabscale;
 300     nouter           = 0;
 301     ninner           = 0;
 302
 303     do
 304     {
 305         #ifdef GMX_THREAD_SHM_FDECOMP
 306         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 307         nn0              = *count;
 308         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 309         *count           = nn1;
 310         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 311         if(nn1>nri) nn1=nri;
 312         #else
 313         nn0 = 0;
 314         nn1 = nri;
 315         #endif
 316
 317         for(n=nn0; (n<nn1); n++)
 318         {
 319             is3              = 3*shift[n];
 320             shX              = shiftvec[is3];
 321             shY              = shiftvec[is3+1];
 322             shZ              = shiftvec[is3+2];
 323             nj0              = jindex[n];
 324             nj1              = jindex[n+1];
 325             ii               = iinr[n];
 326             ii3              = 3*ii;
 327             ix1              = shX + pos[ii3+0];
 328             iy1              = shY + pos[ii3+1];
 329             iz1              = shZ + pos[ii3+2];
 330             iq               = facel*charge[ii];
 331             isai             = invsqrta[ii];
 332             weight_cg1       = wf[ii];
 333             vctot            = 0;
 334             dvdasum          = 0;
 335             fix1             = 0;
 336             fiy1             = 0;
 337             fiz1             = 0;
 338
 339             for(k=nj0; (k<nj1); k++)
 340             {
 341                 jnr              = jjnr[k];
 342                 weight_cg2       = wf[jnr];
 343                 weight_product   = weight_cg1*weight_cg2;
 344                 if (weight_product < ALMOST_ZERO) {
 345                 /* force is zero, skip this molecule */
 346                  continue;
 347                 }
 348                 else if (weight_product >= ALMOST_ONE)
 349                 {
 350                        hybscal = 1.0;
 351                 }
 352                 else
 353                 {
 354                    hybscal = weight_product;
 355                 }
 356                 j3               = 3*jnr;
 357                 jx1              = pos[j3+0];
 358                 jy1              = pos[j3+1];
 359                 jz1              = pos[j3+2];
 360                 dx11             = ix1 - jx1;
 361                 dy11             = iy1 - jy1;
 362                 dz11             = iz1 - jz1;
 363                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 364                 rinv11           = 1.0/sqrt(rsq11);
 365                 isaj             = invsqrta[jnr];
 366                 isaprod          = isai*isaj;
 367                 qq               = iq*charge[jnr];
 368                 vcoul            = qq*rinv11;
 369                 fscal            = vcoul*rinv11;
 370                 qq               = isaprod*(-qq);
 371                 gbscale          = isaprod*gbtabscale;
 372                 dvdaj            = dvda[jnr];
 373                 r                = rsq11*rinv11;
 374                 rt               = r*gbscale;
 375                 n0               = rt;
 376                 eps              = rt-n0;
 377                 eps2             = eps*eps;
 378                 nnn              = 4*n0;
 379                 Y                = GBtab[nnn];
 380                 F                = GBtab[nnn+1];
 381                 Geps             = eps*GBtab[nnn+2];
 382                 Heps2            = eps2*GBtab[nnn+3];
 383                 Fp               = F+Geps+Heps2;
 384                 VV               = Y+eps*Fp;
 385                 FF               = Fp+Geps+2.0*Heps2;
 386                 vgb              = qq*VV;
 387                 fijC             = qq*FF*gbscale;
 388                 dvdatmp          = -0.5*(vgb+fijC*r);
 389                 dvdasum          = dvdasum + dvdatmp;
 390                 dvda[jnr]        = dvdaj+dvdatmp*isaj*isaj;
 391                 vctot            = vctot + vcoul;
 392                 fscal            = -(fijC-fscal)*rinv11;
 393                 fscal *= hybscal;
 394                 tx               = fscal*dx11;
 395                 ty               = fscal*dy11;
 396                 tz               = fscal*dz11;
 397                 fix1             = fix1 + tx;
 398                 fiy1             = fiy1 + ty;
 399                 fiz1             = fiz1 + tz;
 400                 faction[j3+0]    = faction[j3+0] - tx;
 401                 faction[j3+1]    = faction[j3+1] - ty;
 402                 faction[j3+2]    = faction[j3+2] - tz;
 403             }
 404
 405             faction[ii3+0]   = faction[ii3+0] + fix1;
 406             faction[ii3+1]   = faction[ii3+1] + fiy1;
 407             faction[ii3+2]   = faction[ii3+2] + fiz1;
 408             fshift[is3]      = fshift[is3]+fix1;
 409             fshift[is3+1]    = fshift[is3+1]+fiy1;
 410             fshift[is3+2]    = fshift[is3+2]+fiz1;
 411             ggid             = gid[n];
 412             Vc[ggid]         = Vc[ggid] + vctot;
 413             dvda[ii]         = dvda[ii] + dvdasum*isai*isai;
 414             ninner           = ninner + nj1 - nj0;
 415         }
 416
 417         nouter           = nouter + nn1 - nn0;
 418     }
 419     while (nn1<nri);
 420
 421     *outeriter       = nouter;
 422     *inneriter       = ninner;
 423 }
 424
 425