src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel200_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_kernel200_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel200_adress_cg
  33  * Coulomb interaction:     Reaction field
  34  * VdW interaction:         Not calculated
  35  * water optimization:      No
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel200_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         iq;
  80     real         qq,vcoul,vctot;
  81     real         krsq;
  82     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  83     real         jx1,jy1,jz1;
  84     real         dx11,dy11,dz11,rsq11,rinv11;
  85     real         weight_cg1, weight_cg2, weight_product;
  86     real         hybscal;
  87
  88     nri              = *p_nri;
  89     ntype            = *p_ntype;
  90     nthreads         = *p_nthreads;
  91     facel            = *p_facel;
  92     krf              = *p_krf;
  93     crf              = *p_crf;
  94     tabscale         = *p_tabscale;
  95     nouter           = 0;
  96     ninner           = 0;
  97
  98     do
  99     {
 100         #ifdef GMX_THREAD_SHM_FDECOMP
 101         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 102         nn0              = *count;
 103         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 104         *count           = nn1;
 105         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 106         if(nn1>nri) nn1=nri;
 107         #else
 108         nn0 = 0;
 109         nn1 = nri;
 110         #endif
 111
 112         for(n=nn0; (n<nn1); n++)
 113         {
 114             is3              = 3*shift[n];
 115             shX              = shiftvec[is3];
 116             shY              = shiftvec[is3+1];
 117             shZ              = shiftvec[is3+2];
 118             nj0              = jindex[n];
 119             nj1              = jindex[n+1];
 120             ii               = iinr[n];
 121             ii3              = 3*ii;
 122             ix1              = shX + pos[ii3+0];
 123             iy1              = shY + pos[ii3+1];
 124             iz1              = shZ + pos[ii3+2];
 125             iq               = facel*charge[ii];
 126             weight_cg1       = wf[ii];
 127             vctot            = 0;
 128             fix1             = 0;
 129             fiy1             = 0;
 130             fiz1             = 0;
 131
 132             for(k=nj0; (k<nj1); k++)
 133             {
 134                 jnr              = jjnr[k];
 135                 weight_cg2       = wf[jnr];
 136                 weight_product   = weight_cg1*weight_cg2;
 137                 if (weight_product < ALMOST_ZERO) {
 138                        hybscal = 1.0;
 139                 }
 140                 else if (weight_product >= ALMOST_ONE)
 141                 {
 142                   /* force is zero, skip this molecule */
 143                        continue;
 144                 }
 145                 else
 146                 {
 147                    hybscal = 1.0 - weight_product;
 148                 }
 149                 j3               = 3*jnr;
 150                 jx1              = pos[j3+0];
 151                 jy1              = pos[j3+1];
 152                 jz1              = pos[j3+2];
 153                 dx11             = ix1 - jx1;
 154                 dy11             = iy1 - jy1;
 155                 dz11             = iz1 - jz1;
 156                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 157                 rinv11           = 1.0/sqrt(rsq11);
 158                 qq               = iq*charge[jnr];
 159                 rinvsq           = rinv11*rinv11;
 160                 krsq             = krf*rsq11;
 161                 vcoul            = qq*(rinv11+krsq-crf);
 162                 vctot            = vctot+vcoul;
 163                 fscal            = (qq*(rinv11-2.0*krsq))*rinvsq;
 164                 fscal *= hybscal;
 165                 tx               = fscal*dx11;
 166                 ty               = fscal*dy11;
 167                 tz               = fscal*dz11;
 168                 fix1             = fix1 + tx;
 169                 fiy1             = fiy1 + ty;
 170                 fiz1             = fiz1 + tz;
 171                 faction[j3+0]    = faction[j3+0] - tx;
 172                 faction[j3+1]    = faction[j3+1] - ty;
 173                 faction[j3+2]    = faction[j3+2] - tz;
 174             }
 175
 176             faction[ii3+0]   = faction[ii3+0] + fix1;
 177             faction[ii3+1]   = faction[ii3+1] + fiy1;
 178             faction[ii3+2]   = faction[ii3+2] + fiz1;
 179             fshift[is3]      = fshift[is3]+fix1;
 180             fshift[is3+1]    = fshift[is3+1]+fiy1;
 181             fshift[is3+2]    = fshift[is3+2]+fiz1;
 182             ggid             = gid[n];
 183             Vc[ggid]         = Vc[ggid] + vctot;
 184             ninner           = ninner + nj1 - nj0;
 185         }
 186
 187         nouter           = nouter + nn1 - nn0;
 188     }
 189     while (nn1<nri);
 190
 191     *outeriter       = nouter;
 192     *inneriter       = ninner;
 193 }
 194
 195
 196
 197
 198
 199 /*
 200  * Gromacs nonbonded kernel nb_kernel200_adress_ex
 201  * Coulomb interaction:     Reaction field
 202  * VdW interaction:         Not calculated
 203  * water optimization:      No
 204  * Calculate forces:        yes
 205  */
 206 void nb_kernel200_adress_ex(
 207                     int *           p_nri,
 208                     int *           iinr,
 209                     int *           jindex,
 210                     int *           jjnr,
 211                     int *           shift,
 212                     real *         shiftvec,
 213                     real *         fshift,
 214                     int *           gid,
 215                     real *         pos,
 216                     real *         faction,
 217                     real *         charge,
 218                     real *         p_facel,
 219                     real *         p_krf,
 220                     real *         p_crf,
 221                     real *         Vc,
 222                     int *           type,
 223                     int *           p_ntype,
 224                     real *         vdwparam,
 225                     real *         Vvdw,
 226                     real *         p_tabscale,
 227                     real *         VFtab,
 228                     real *         invsqrta,
 229                     real *         dvda,
 230                     real *         p_gbtabscale,
 231                     real *         GBtab,
 232                     int *           p_nthreads,
 233                     int *           count,
 234                     void *          mtx,
 235                     int *           outeriter,
 236                     int *           inneriter,
 237                     real           force_cap,
 238                     real *         wf)
 239 {
 240     int           nri,ntype,nthreads;
 241     real         facel,krf,crf,tabscale,gbtabscale;
 242     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 243     int           nn0,nn1,nouter,ninner;
 244     real         shX,shY,shZ;
 245     real         fscal,tx,ty,tz;
 246     real         rinvsq;
 247     real         iq;
 248     real         qq,vcoul,vctot;
 249     real         krsq;
 250     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 251     real         jx1,jy1,jz1;
 252     real         dx11,dy11,dz11,rsq11,rinv11;
 253     real         weight_cg1, weight_cg2, weight_product;
 254     real         hybscal;
 255
 256     nri              = *p_nri;
 257     ntype            = *p_ntype;
 258     nthreads         = *p_nthreads;
 259     facel            = *p_facel;
 260     krf              = *p_krf;
 261     crf              = *p_crf;
 262     tabscale         = *p_tabscale;
 263     nouter           = 0;
 264     ninner           = 0;
 265
 266     do
 267     {
 268         #ifdef GMX_THREAD_SHM_FDECOMP
 269         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 270         nn0              = *count;
 271         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 272         *count           = nn1;
 273         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 274         if(nn1>nri) nn1=nri;
 275         #else
 276         nn0 = 0;
 277         nn1 = nri;
 278         #endif
 279
 280         for(n=nn0; (n<nn1); n++)
 281         {
 282             is3              = 3*shift[n];
 283             shX              = shiftvec[is3];
 284             shY              = shiftvec[is3+1];
 285             shZ              = shiftvec[is3+2];
 286             nj0              = jindex[n];
 287             nj1              = jindex[n+1];
 288             ii               = iinr[n];
 289             ii3              = 3*ii;
 290             ix1              = shX + pos[ii3+0];
 291             iy1              = shY + pos[ii3+1];
 292             iz1              = shZ + pos[ii3+2];
 293             iq               = facel*charge[ii];
 294             weight_cg1       = wf[ii];
 295             vctot            = 0;
 296             fix1             = 0;
 297             fiy1             = 0;
 298             fiz1             = 0;
 299
 300             for(k=nj0; (k<nj1); k++)
 301             {
 302                 jnr              = jjnr[k];
 303                 weight_cg2       = wf[jnr];
 304                 weight_product   = weight_cg1*weight_cg2;
 305                 if (weight_product < ALMOST_ZERO) {
 306                 /* force is zero, skip this molecule */
 307                  continue;
 308                 }
 309                 else if (weight_product >= ALMOST_ONE)
 310                 {
 311                        hybscal = 1.0;
 312                 }
 313                 else
 314                 {
 315                    hybscal = weight_product;
 316                 }
 317                 j3               = 3*jnr;
 318                 jx1              = pos[j3+0];
 319                 jy1              = pos[j3+1];
 320                 jz1              = pos[j3+2];
 321                 dx11             = ix1 - jx1;
 322                 dy11             = iy1 - jy1;
 323                 dz11             = iz1 - jz1;
 324                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 325                 rinv11           = 1.0/sqrt(rsq11);
 326                 qq               = iq*charge[jnr];
 327                 rinvsq           = rinv11*rinv11;
 328                 krsq             = krf*rsq11;
 329                 vcoul            = qq*(rinv11+krsq-crf);
 330                 vctot            = vctot+vcoul;
 331                 fscal            = (qq*(rinv11-2.0*krsq))*rinvsq;
 332                 fscal *= hybscal;
 333                 tx               = fscal*dx11;
 334                 ty               = fscal*dy11;
 335                 tz               = fscal*dz11;
 336                 fix1             = fix1 + tx;
 337                 fiy1             = fiy1 + ty;
 338                 fiz1             = fiz1 + tz;
 339                 faction[j3+0]    = faction[j3+0] - tx;
 340                 faction[j3+1]    = faction[j3+1] - ty;
 341                 faction[j3+2]    = faction[j3+2] - tz;
 342             }
 343
 344             faction[ii3+0]   = faction[ii3+0] + fix1;
 345             faction[ii3+1]   = faction[ii3+1] + fiy1;
 346             faction[ii3+2]   = faction[ii3+2] + fiz1;
 347             fshift[is3]      = fshift[is3]+fix1;
 348             fshift[is3+1]    = fshift[is3+1]+fiy1;
 349             fshift[is3+2]    = fshift[is3+2]+fiz1;
 350             ggid             = gid[n];
 351             Vc[ggid]         = Vc[ggid] + vctot;
 352             ninner           = ninner + nj1 - nj0;
 353         }
 354
 355         nouter           = nouter + nn1 - nn0;
 356     }
 357     while (nn1<nri);
 358
 359     *outeriter       = nouter;
 360     *inneriter       = ninner;
 361 }
 362
 363