src/gmxlib/nonbonded/nb_generic_adress.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43
  44 #include "types/simple.h"
  45 #include "vec.h"
  46 #include "typedefs.h"
  47 #include "nb_generic_adress.h"
  48
  49 #define ALMOST_ZERO 1e-30
  50 #define ALMOST_ONE 1-(1e-30)
  51 void
  52 gmx_nb_generic_adress_kernel(t_nblist *           nlist,
  53                              t_forcerec *         fr,
  54                              t_mdatoms *          mdatoms,
  55                              real *               x,
  56                              real *               f,
  57                              real *               fshift,
  58                              real *               Vc,
  59                              real *               Vvdw,
  60                              real                 tabscale,
  61                              real *               VFtab,
  62                              int *                outeriter,
  63                              int *                inneriter,
  64                              gmx_bool                bCG)
  65 {
  66     int           nri,ntype,table_nelements,ielec,ivdw;
  67     real          facel,gbtabscale;
  68     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid,nnn,n0;
  69     real          shX,shY,shZ;
  70     real          fscal,tx,ty,tz;
  71     real          rinvsq;
  72     real          iq;
  73     real          qq,vcoul,krsq,vctot;
  74     int           nti,nvdwparam;
  75     int           tj;
  76     real          rt,r,eps,eps2,Y,F,Geps,Heps2,VV,FF,Fp,fijD,fijR;
  77     real          rinvsix;
  78     real          Vvdwtot;
  79     real          Vvdw_rep,Vvdw_disp;
  80     real          ix,iy,iz,fix,fiy,fiz;
  81     real          jx,jy,jz;
  82     real          dx,dy,dz,rsq,rinv;
  83     real          c6,c12,cexp1,cexp2,br;
  84     real *        charge;
  85     real *        shiftvec;
  86     real *        vdwparam;
  87     int *         shift;
  88     int *         type;
  89
  90     real *     wf;
  91     real       weight_cg1;
  92     real       weight_cg2;
  93     real       weight_product;
  94     real       hybscal; /* the multiplicator to the force for hybrid interactions*/
  95     gmx_bool   bHybrid; /*Are we in the hybrid zone ?*/
  96     real       force_cap;
  97
  98     wf                  = mdatoms->wf;
  99
 100     force_cap = fr->adress_ex_forcecap;
 101
 102     ielec               = nlist->ielec;
 103     ivdw                = nlist->ivdw;
 104
 105     /* avoid compiler warnings for cases that cannot happen */
 106     nnn                 = 0;
 107     vcoul               = 0.0;
 108     eps                 = 0.0;
 109     eps2                = 0.0;
 110
 111     /* 3 VdW parameters for buckingham, otherwise 2 */
 112     nvdwparam           = (nlist->ivdw==2) ? 3 : 2;
 113     table_nelements     = (ielec==3) ? 4 : 0;
 114     table_nelements    += (ivdw==3) ? 8 : 0;
 115
 116     charge              = mdatoms->chargeA;
 117     type                = mdatoms->typeA;
 118     facel               = fr->epsfac;
 119     shiftvec            = fr->shift_vec[0];
 120     vdwparam            = fr->nbfp;
 121     ntype               = fr->ntype;
 122
 123
 124
 125
 126    for(n=0; (n<nlist->nri); n++)
 127     {
 128         is3              = 3*nlist->shift[n];
 129         shX              = shiftvec[is3];
 130         shY              = shiftvec[is3+1];
 131         shZ              = shiftvec[is3+2];
 132         nj0              = nlist->jindex[n];
 133         nj1              = nlist->jindex[n+1];
 134         ii               = nlist->iinr[n];
 135         ii3              = 3*ii;
 136         ix               = shX + x[ii3+0];
 137         iy               = shY + x[ii3+1];
 138         iz               = shZ + x[ii3+2];
 139         iq               = facel*charge[ii];
 140         nti              = nvdwparam*ntype*type[ii];
 141         vctot            = 0;
 142         Vvdwtot          = 0;
 143         fix              = 0;
 144         fiy              = 0;
 145         fiz              = 0;
 146
 147         weight_cg1       = wf[ii];
 148
 149         /* TODO: why does this line her not speed up things ?
 150          * if ((!bCG) && weight_cg1 < ALMOST_ZERO) continue;
 151          */
 152         for(k=nj0; (k<nj1); k++)
 153         {
 154             jnr              = nlist->jjnr[k];
 155             weight_cg2       = wf[jnr];
 156
 157             weight_product   = weight_cg1*weight_cg2;
 158
 159             if (weight_product < ALMOST_ZERO)
 160             {
 161                 /* if it's a explicit loop, skip this atom */
 162                 if (!bCG)
 163                 {
 164                     continue;
 165                 }
 166                 else /* if it's a coarse grained loop, include this atom */
 167                 {
 168                     bHybrid = FALSE;
 169                     hybscal = 1.0;
 170                 }
 171             }
 172             else if (weight_product >= ALMOST_ONE)
 173             {
 174
 175                 /* if it's a explicit loop, include this atom */
 176                 if(!bCG)
 177                 {
 178                     bHybrid = FALSE;
 179                     hybscal = 1.0;
 180                 }
 181                 else  /* if it's a coarse grained loop, skip this atom */
 182                 {
 183                     continue;
 184                 }
 185             }
 186             /* both have double identity, get hybrid scaling factor */
 187             else
 188             {
 189                 bHybrid = TRUE;
 190                 hybscal = weight_product;
 191
 192                 if(bCG)
 193                 {
 194                     hybscal = 1.0 - hybscal;
 195                 }
 196             }
 197
 198
 199             j3               = 3*jnr;
 200             jx               = x[j3+0];
 201             jy               = x[j3+1];
 202             jz               = x[j3+2];
 203             dx               = ix - jx;
 204             dy               = iy - jy;
 205             dz               = iz - jz;
 206             rsq              = dx*dx+dy*dy+dz*dz;
 207             rinv             = gmx_invsqrt(rsq);
 208             rinvsq           = rinv*rinv;
 209
 210
 211             fscal            = 0;
 212
 213             if(ielec==3 || ivdw==3)
 214             {
 215                 r                = rsq*rinv;
 216                 rt               = r*tabscale;
 217                 n0               = rt;
 218                 eps              = rt-n0;
 219                 eps2             = eps*eps;
 220                 nnn              = table_nelements*n0;
 221             }
 222
 223             /* Coulomb interaction. ielec==0 means no interaction */
 224             if(ielec>0)
 225             {
 226                 qq               = iq*charge[jnr];
 227
 228                 switch(ielec)
 229                 {
 230                     case 1:
 231                         /* Vanilla cutoff coulomb */
 232                         vcoul            = qq*rinv;
 233                         fscal            = vcoul*rinvsq;
 234                         break;
 235
 236                     case 2:
 237                         /* Reaction-field */
 238                         krsq             = fr->k_rf*rsq;
 239                         vcoul            = qq*(rinv+krsq-fr->c_rf);
 240                         fscal            = qq*(rinv-2.0*krsq)*rinvsq;
 241                         break;
 242
 243                     case 3:
 244                         /* Tabulated coulomb */
 245                         Y                = VFtab[nnn];
 246                         F                = VFtab[nnn+1];
 247                         Geps             = eps*VFtab[nnn+2];
 248                         Heps2            = eps2*VFtab[nnn+3];
 249                         nnn             += 4;
 250                         Fp               = F+Geps+Heps2;
 251                         VV               = Y+eps*Fp;
 252                         FF               = Fp+Geps+2.0*Heps2;
 253                         vcoul            = qq*VV;
 254                         fscal            = -qq*FF*tabscale*rinv;
 255                         break;
 256
 257                     case 4:
 258                         /* GB */
 259                         gmx_fatal(FARGS,"Death & horror! GB generic interaction not implemented.\n");
 260                         break;
 261
 262                     default:
 263                         gmx_fatal(FARGS,"Death & horror! No generic coulomb interaction for ielec=%d.\n",ielec);
 264                         break;
 265                 }
 266                 vctot            = vctot+vcoul;
 267             } /* End of coulomb interactions */
 268
 269             /* VdW interaction. ivdw==0 means no interaction */
 270             if(ivdw>0)
 271             {
 272                 tj               = nti+nvdwparam*type[jnr];
 273
 274                 switch(ivdw)
 275                 {
 276                                         case 1:
 277                                                 /* Vanilla Lennard-Jones cutoff */
 278                                                 c6               = vdwparam[tj];
 279                                                 c12              = vdwparam[tj+1];
 280
 281                                                 rinvsix          = rinvsq*rinvsq*rinvsq;
 282                                                 Vvdw_disp        = c6*rinvsix;
 283                                                 Vvdw_rep         = c12*rinvsix*rinvsix;
 284                                                 fscal           += (12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
 285                                                 Vvdwtot          = Vvdwtot+Vvdw_rep-Vvdw_disp;
 286                                                 break;
 287
 288                                         case 2:
 289                                                 /* Buckingham */
 290                                                 c6               = vdwparam[tj];
 291                                                 cexp1            = vdwparam[tj+1];
 292                                                 cexp2            = vdwparam[tj+2];
 293
 294                                                 rinvsix          = rinvsq*rinvsq*rinvsq;
 295                                                 Vvdw_disp        = c6*rinvsix;
 296                                                 br               = cexp2*rsq*rinv;
 297                                                 Vvdw_rep         = cexp1*exp(-br);
 298                                                 fscal           += (br*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
 299                                                 Vvdwtot          = Vvdwtot+Vvdw_rep-Vvdw_disp;
 300                                                 break;
 301
 302                                         case 3:
 303                                                 /* Tabulated VdW */
 304                                                 c6               = vdwparam[tj];
 305                                                 c12              = vdwparam[tj+1];
 306
 307                                                 Y                = VFtab[nnn];
 308                                                 F                = VFtab[nnn+1];
 309                                                 Geps             = eps*VFtab[nnn+2];
 310                                                 Heps2            = eps2*VFtab[nnn+3];
 311                                                 Fp               = F+Geps+Heps2;
 312                                                 VV               = Y+eps*Fp;
 313                                                 FF               = Fp+Geps+2.0*Heps2;
 314                                                 Vvdw_disp        = c6*VV;
 315                                                 fijD             = c6*FF;
 316                                                 nnn             += 4;
 317                                                 Y                = VFtab[nnn];
 318                                                 F                = VFtab[nnn+1];
 319                                                 Geps             = eps*VFtab[nnn+2];
 320                                                 Heps2            = eps2*VFtab[nnn+3];
 321                                                 Fp               = F+Geps+Heps2;
 322                                                 VV               = Y+eps*Fp;
 323                                                 FF               = Fp+Geps+2.0*Heps2;
 324                                                 Vvdw_rep         = c12*VV;
 325                                                 fijR             = c12*FF;
 326                                                 fscal           += -(fijD+fijR)*tabscale*rinv;
 327                                                 Vvdwtot          = Vvdwtot + Vvdw_disp + Vvdw_rep;
 328                                                 if(!bCG && force_cap>0 && (fabs(fscal)> force_cap))
 329                                                 {
 330                                                      fscal=force_cap*fscal/fabs(fscal);
 331                                                 }
 332                                                 break;
 333
 334                                         default:
 335                                                 gmx_fatal(FARGS,"Death & horror! No generic VdW interaction for ivdw=%d.\n",ivdw);
 336                                                 break;
 337                                 }
 338                         } /* end VdW interactions */
 339
 340              /* force weight is one anyway */
 341                     if (bHybrid)
 342                     {
 343                         fscal *= hybscal;
 344                     }
 345
 346             tx               = fscal*dx;
 347             ty               = fscal*dy;
 348             tz               = fscal*dz;
 349             fix              = fix + tx;
 350             fiy              = fiy + ty;
 351             fiz              = fiz + tz;
 352             f[j3+0]          = f[j3+0] - tx;
 353             f[j3+1]          = f[j3+1] - ty;
 354             f[j3+2]          = f[j3+2] - tz;
 355         }
 356
 357         f[ii3+0]         = f[ii3+0] + fix;
 358         f[ii3+1]         = f[ii3+1] + fiy;
 359         f[ii3+2]         = f[ii3+2] + fiz;
 360         fshift[is3]      = fshift[is3]+fix;
 361         fshift[is3+1]    = fshift[is3+1]+fiy;
 362         fshift[is3+2]    = fshift[is3+2]+fiz;
 363         ggid             = nlist->gid[n];
 364         Vc[ggid]         = Vc[ggid] + vctot;
 365         Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 366     }
 367
 368     *outeriter       = nlist->nri;
 369     *inneriter       = nlist->jindex[n];
 370 }
 371