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