src/gmxlib/nonbonded/nb_free_energy.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41
  42 #include "vec.h"
  43 #include "typedefs.h"
  44
  45 void
  46 gmx_nb_free_energy_kernel(int                  icoul,
  47                           int                  ivdw,
  48                           int                  nri,
  49                           int *                iinr,
  50                           int *                jindex,
  51                           int *                jjnr,
  52                           int *                shift,
  53                           real *               shiftvec,
  54                           real *               fshift,
  55                           int *                gid,
  56                           real *               x,
  57                           real *               f,
  58                           real *               chargeA,
  59                           real *               chargeB,
  60                           real                 facel,
  61                           real                 krf,
  62                           real                 crf,
  63                           real                 ewc,
  64                           real *               Vc,
  65                           int *                typeA,
  66                           int *                typeB,
  67                           int                  ntype,
  68                           real *               nbfp,
  69                           real *               Vvdw,
  70                           real                 tabscale,
  71                           real *               VFtab,
  72                           real                 lambda,
  73                           real *               dvdlambda,
  74                           real                 alpha,
  75                           int                  lam_power,
  76                           real                 sigma6_def,
  77                           real                 sigma6_min,
  78                           gmx_bool                 bDoForces,
  79                           int *                outeriter,
  80                           int *                inneriter)
  81 {
  82     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
  83     real          shX,shY,shZ;
  84     real          Fscal,FscalA,FscalB,tx,ty,tz;
  85     real          VcoulA,VcoulB,VvdwA,VvdwB;
  86     real          rinv6,r,rt;
  87     real          iqA,iqB;
  88     real          qqA,qqB,vcoul,vctot,krsq;
  89     int           ntiA,ntiB;
  90     int           tjA,tjB;
  91     real          rinvsix;
  92     real          Vvdw6,Vvdwtot;
  93     real          Vvdw12;
  94     real          ix,iy,iz,fix,fiy,fiz;
  95     real          dx,dy,dz,rsq,r4,r6,rinv;
  96     real          c6A,c12A,c6B,c12B;
  97     real          dvdl,L1,alfA,alfB,dalfA,dalfB;
  98     real          sigma6a,sigma6b;
  99     real          rA,rinvA,rinv4A,rB,rinvB,rinv4B;
 100     int           do_coultab,do_vdwtab,do_tab,tab_elemsize;
 101     int           n0,n1,nnn;
 102     real          Y,F,G,H,Fp,Geps,Heps2,eps,eps2,VV,FF;
 103     double        isp=0.564189583547756;
 104
 105
 106     /* fix compiler warnings */
 107     nj1 = 0;
 108     n1  = 0;
 109     eps = 0;
 110     eps2 = 0;
 111
 112     dvdl = 0;
 113     L1   = 1.0 - lambda;
 114
 115     alfA  = alpha*(lam_power==2 ? lambda*lambda : lambda);
 116     alfB  = alpha*(lam_power==2 ? L1*L1 : L1);
 117     dalfA = alpha*lam_power/6.0*(lam_power==2 ? lambda : 1);
 118     dalfB = alpha*lam_power/6.0*(lam_power==2 ? L1 : 1);
 119
 120     /* Ewald table is special (icoul==5) */
 121
 122     do_coultab = (icoul==3);
 123     do_vdwtab  = (ivdw==3);
 124
 125     do_tab = do_coultab || do_vdwtab;
 126
 127     /* we always use the combined table here */
 128     tab_elemsize = 12;
 129
 130     for(n=0; (n<nri); n++)
 131     {
 132         is3              = 3*shift[n];
 133         shX              = shiftvec[is3];
 134         shY              = shiftvec[is3+1];
 135         shZ              = shiftvec[is3+2];
 136         nj0              = jindex[n];
 137         nj1              = jindex[n+1];
 138         ii               = iinr[n];
 139         ii3              = 3*ii;
 140         ix               = shX + x[ii3+0];
 141         iy               = shY + x[ii3+1];
 142         iz               = shZ + x[ii3+2];
 143         iqA              = facel*chargeA[ii];
 144         iqB              = facel*chargeB[ii];
 145         ntiA             = 2*ntype*typeA[ii];
 146         ntiB             = 2*ntype*typeB[ii];
 147         vctot            = 0;
 148         Vvdwtot          = 0;
 149         fix              = 0;
 150         fiy              = 0;
 151         fiz              = 0;
 152
 153         for(k=nj0; (k<nj1); k++)
 154         {
 155             jnr              = jjnr[k];
 156             j3               = 3*jnr;
 157             dx               = ix - x[j3];
 158             dy               = iy - x[j3+1];
 159             dz               = iz - x[j3+2];
 160             rsq              = dx*dx+dy*dy+dz*dz;
 161             rinv             = gmx_invsqrt(rsq);
 162             r                = rsq*rinv;
 163             tjA              = ntiA+2*typeA[jnr];
 164             tjB              = ntiB+2*typeB[jnr];
 165             c6A              = nbfp[tjA];
 166             c6B              = nbfp[tjB];
 167             c12A             = nbfp[tjA+1];
 168             c12B             = nbfp[tjB+1];
 169             qqA              = iqA*chargeA[jnr];
 170             qqB              = iqB*chargeB[jnr];
 171
 172             if((c6A > 0) && (c12A > 0))
 173             {
 174                 sigma6a      = c12A/c6A;
 175
 176                 if (sigma6a < sigma6_min)
 177                 {
 178                     sigma6a  = sigma6_min;
 179                 }
 180             }
 181             else
 182             {
 183                 sigma6a      = sigma6_def;
 184             }
 185             if((c6B > 0) && (c12B > 0))
 186             {
 187                 sigma6b      = c12B/c6B;
 188
 189                 if (sigma6b < sigma6_min)
 190                 {
 191                     sigma6b  = sigma6_min;
 192                 }
 193             }
 194             else
 195             {
 196                 sigma6b      = sigma6_def;
 197             }
 198
 199             r4               = rsq*rsq;
 200             r6               = r4*rsq;
 201
 202             FscalA           = 0;
 203             VcoulA           = 0;
 204             VvdwA            = 0;
 205             rinv4A           = 0;
 206
 207             /* Only spend time on A state if it is non-zero */
 208             if( (qqA != 0) || (c6A != 0) || (c12A != 0) )
 209             {
 210                 rA             = pow(alfA*sigma6a+r6,1.0/6.0);
 211                 rinvA          = 1.0/rA;
 212                 rinv4A         = rinvA*rinvA;
 213                 rinv4A         = rinv4A*rinv4A;
 214
 215
 216                 if(do_tab)
 217                 {
 218                     rt         = rA*tabscale;
 219                     n0         = rt;
 220                     eps        = rt-n0;
 221                     eps2       = eps*eps;
 222                     n1         = tab_elemsize*n0;
 223                 }
 224
 225                 if(icoul==1 || icoul==5)
 226                 {
 227                     /* simple cutoff */
 228                     VcoulA     = qqA*rinvA;
 229                     FscalA     = VcoulA*rinvA*rinvA;
 230                 }
 231                 else if(icoul==2)
 232                 {
 233                     /* reaction-field */
 234                     krsq       = krf*rA*rA;
 235                     VcoulA     = qqA*(rinvA+krsq-crf);
 236                     FscalA     = qqA*(rinvA-2.0*krsq)*rinvA*rinvA;
 237                 }
 238                 else if(icoul==3)
 239                 {
 240                     /* non-Ewald tabulated coulomb */
 241                     nnn        = n1;
 242                     Y          = VFtab[nnn];
 243                     F          = VFtab[nnn+1];
 244                     Geps       = eps*VFtab[nnn+2];
 245                     Heps2      = eps2*VFtab[nnn+3];
 246                     Fp         = F+Geps+Heps2;
 247                     VV         = Y+eps*Fp;
 248                     FF         = Fp+Geps+2.0*Heps2;
 249                     VcoulA     = qqA*VV;
 250                     FscalA     = -qqA*tabscale*FF*rinvA;
 251                 }
 252
 253                 if(ivdw==1)
 254                 {
 255                     /* cutoff LJ */
 256                     rinv6            = rinvA*rinvA*rinv4A;
 257                     Vvdw6            = c6A*rinv6;
 258                     Vvdw12           = c12A*rinv6*rinv6;
 259                     VvdwA            = Vvdw12-Vvdw6;
 260                     FscalA          += (12.0*Vvdw12-6.0*Vvdw6)*rinvA*rinvA;
 261                 }
 262                 else if(ivdw==3)
 263                 {
 264                     /* Table LJ */
 265                     nnn = n1+4;
 266
 267                     /* dispersion */
 268                     Y          = VFtab[nnn];
 269                     F          = VFtab[nnn+1];
 270                     Geps       = eps*VFtab[nnn+2];
 271                     Heps2      = eps2*VFtab[nnn+3];
 272                     Fp         = F+Geps+Heps2;
 273                     VV         = Y+eps*Fp;
 274                     FF         = Fp+Geps+2.0*Heps2;
 275                     VvdwA     += c6A*VV;
 276                     FscalA    -= c6A*tabscale*FF*rinvA;
 277
 278                     /* repulsion */
 279                     Y          = VFtab[nnn+4];
 280                     F          = VFtab[nnn+5];
 281                     Geps       = eps*VFtab[nnn+6];
 282                     Heps2      = eps2*VFtab[nnn+7];
 283                     Fp         = F+Geps+Heps2;
 284                     VV         = Y+eps*Fp;
 285                     FF         = Fp+Geps+2.0*Heps2;
 286                     VvdwA     += c12A*VV;
 287                     FscalA    -= c12A*tabscale*FF*rinvA;
 288                 }
 289                 /* Buckingham vdw free energy not supported */
 290             }
 291
 292             FscalB           = 0;
 293             VcoulB           = 0;
 294             VvdwB            = 0;
 295             rinv4B           = 0;
 296
 297             /* Only spend time on B state if it is non-zero */
 298             if( (qqB != 0) || (c6B != 0) || (c12B != 0) )
 299             {
 300                 rB             = pow(alfB*sigma6b+r6,1.0/6.0);
 301                 rinvB          = 1.0/rB;
 302                 rinv4B         = rinvB*rinvB;
 303                 rinv4B         = rinv4B*rinv4B;
 304
 305
 306                 if(do_tab)
 307                 {
 308                     rt         = rB*tabscale;
 309                     n0         = rt;
 310                     eps        = rt-n0;
 311                     eps2       = eps*eps;
 312                     n1         = tab_elemsize*n0;
 313                 }
 314
 315                 if(icoul==1 || icoul==5)
 316                 {
 317                     /* simple cutoff */
 318                     VcoulB     = qqB*rinvB;
 319                     FscalB     = VcoulB*rinvB*rinvB;
 320                 }
 321                 else if(icoul==2)
 322                 {
 323                     /* reaction-field */
 324                     krsq       = krf*rB*rB;
 325                     VcoulB     = qqB*(rinvB+krsq-crf);
 326                     FscalB     = qqB*(rinvB-2.0*krsq)*rinvB*rinvB;
 327                 }
 328                 else if(icoul==3)
 329                 {
 330                     /* non-Ewald tabulated coulomb */
 331                     nnn        = n1;
 332                     Y          = VFtab[nnn];
 333                     F          = VFtab[nnn+1];
 334                     Geps       = eps*VFtab[nnn+2];
 335                     Heps2      = eps2*VFtab[nnn+3];
 336                     Fp         = F+Geps+Heps2;
 337                     VV         = Y+eps*Fp;
 338                     FF         = Fp+Geps+2.0*Heps2;
 339                     VcoulB     = qqB*VV;
 340                     FscalB     = -qqB*tabscale*FF*rinvB;
 341                 }
 342
 343                 if(ivdw==1)
 344                 {
 345                     /* cutoff LJ */
 346                     rinv6            = rinvB*rinvB*rinv4B;
 347                     Vvdw6            = c6B*rinv6;
 348                     Vvdw12           = c12B*rinv6*rinv6;
 349                     VvdwB            = Vvdw12-Vvdw6;
 350                     FscalB          += (12.0*Vvdw12-6.0*Vvdw6)*rinvB*rinvB;
 351                 }
 352                 else if(ivdw==3)
 353                 {
 354                     /* Table LJ */
 355                     nnn = n1+4;
 356
 357                     /* dispersion */
 358                     Y          = VFtab[nnn];
 359                     F          = VFtab[nnn+1];
 360                     Geps       = eps*VFtab[nnn+2];
 361                     Heps2      = eps2*VFtab[nnn+3];
 362                     Fp         = F+Geps+Heps2;
 363                     VV         = Y+eps*Fp;
 364                     FF         = Fp+Geps+2.0*Heps2;
 365                     VvdwB     += c6B*VV;
 366                     FscalB    -= c6B*tabscale*FF*rinvB;
 367
 368                     /* repulsion */
 369                     Y          = VFtab[nnn+4];
 370                     F          = VFtab[nnn+5];
 371                     Geps       = eps*VFtab[nnn+6];
 372                     Heps2      = eps2*VFtab[nnn+7];
 373                     Fp         = F+Geps+Heps2;
 374                     VV         = Y+eps*Fp;
 375                     FF         = Fp+Geps+2.0*Heps2;
 376                     VvdwB     += c12B*VV;
 377                     FscalB    -= c12B*tabscale*FF*rinvB;
 378                 }
 379                 /* Buckingham vdw free energy not supported */
 380             }
 381
 382             Fscal = 0;
 383
 384             if(icoul==5)
 385             {
 386                 /* Soft-core Ewald interactions are special:
 387                  * For the direct space interactions we effectively want the
 388                  * normal coulomb interaction (added above when icoul==5),
 389                  * but need to subtract the part added in reciprocal space.
 390                  */
 391                 if (r != 0)
 392                 {
 393                     VV    = gmx_erf(ewc*r)*rinv;
 394                     FF    = rinv*rinv*(VV - 2.0*ewc*isp*exp(-ewc*ewc*rsq));
 395                 }
 396                 else
 397                 {
 398                     VV    = ewc*2.0/sqrt(M_PI);
 399                     FF    = 0;
 400                 }
 401                 vctot  -= (lambda*qqB + L1*qqA)*VV;
 402                 Fscal  -= (lambda*qqB + L1*qqA)*FF;
 403                 dvdl   -= (qqB - qqA)*VV;
 404             }
 405
 406             /* Assemble A and B states */
 407             vctot         += lambda*VcoulB + L1*VcoulA;
 408             Vvdwtot       += lambda*VvdwB  + L1*VvdwA;
 409
 410             Fscal         += (L1*FscalA*rinv4A + lambda*FscalB*rinv4B)*r4;
 411             dvdl          += (VcoulB + VvdwB) - (VcoulA + VvdwA);
 412             dvdl          += lambda*dalfB*FscalB*sigma6b*rinv4B
 413                                - L1*dalfA*FscalA*sigma6a*rinv4A;
 414
 415             if (bDoForces)
 416             {
 417                 tx         = Fscal*dx;
 418                 ty         = Fscal*dy;
 419                 tz         = Fscal*dz;
 420                 fix        = fix + tx;
 421                 fiy        = fiy + ty;
 422                 fiz        = fiz + tz;
 423                 f[j3]      = f[j3]   - tx;
 424                 f[j3+1]    = f[j3+1] - ty;
 425                 f[j3+2]    = f[j3+2] - tz;
 426             }
 427         }
 428
 429         if (bDoForces)
 430         {
 431             f[ii3]         = f[ii3]        + fix;
 432             f[ii3+1]       = f[ii3+1]      + fiy;
 433             f[ii3+2]       = f[ii3+2]      + fiz;
 434             fshift[is3]    = fshift[is3]   + fix;
 435             fshift[is3+1]  = fshift[is3+1] + fiy;
 436             fshift[is3+2]  = fshift[is3+2] + fiz;
 437         }
 438         ggid               = gid[n];
 439         Vc[ggid]           = Vc[ggid] + vctot;
 440         Vvdw[ggid]         = Vvdw[ggid] + Vvdwtot;
 441     }
 442
 443     *dvdlambda      += dvdl;
 444     *outeriter       = nri;
 445     *inneriter       = nj1;
 446 }