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