src/gromacs/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 *               Vv,
  70                           real                 tabscale,
  71                           real *               VFtab,
  72                           real                 lambda_coul,
  73                           real                 lambda_vdw,
  74                           real *               dvdl,
  75                           real                 alpha_coul,
  76                           real                 alpha_vdw,
  77                           int                  lam_power,
  78                           real                 sc_r_power,
  79                           real                 sigma6_def,
  80                           real                 sigma6_min,
  81                           gmx_bool             bDoForces,
  82                           int *                outeriter,
  83                           int *                inneriter)
  84 {
  85
  86 #define  STATE_A  0
  87 #define  STATE_B  1
  88 #define  NSTATES  2
  89     int           i,j,n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
  90     real          shX,shY,shZ;
  91     real          Fscal,FscalC[NSTATES],FscalV[NSTATES],tx,ty,tz;
  92     real          Vcoul[NSTATES],Vvdw[NSTATES];
  93     real          rinv6,r,rt,rtC,rtV;
  94     real          iqA,iqB;
  95     real          qq[NSTATES],vctot,krsq;
  96     int           ntiA,ntiB,tj[NSTATES];
  97     real          Vvdw6, Vvdw12,vvtot;
  98     real          ix,iy,iz,fix,fiy,fiz;
  99     real          dx,dy,dz,rsq,rinv;
 100     real          c6[NSTATES],c12[NSTATES];
 101     real          LFC[NSTATES],LFV[NSTATES],DLF[NSTATES];
 102     double        dvdl_coul,dvdl_vdw;
 103     real          lfac_coul[NSTATES],dlfac_coul[NSTATES],lfac_vdw[NSTATES],dlfac_vdw[NSTATES];
 104     real          sigma6[NSTATES],alpha_vdw_eff,alpha_coul_eff,sigma2_def,sigma2_min;
 105     real          rp,rpm2,rC,rV,rinvC,rpinvC,rinvV,rpinvV;
 106     real          sigma2[NSTATES],sigma_pow[NSTATES],sigma_powm2[NSTATES],rs,rs2;
 107     int           do_coultab,do_vdwtab,do_tab,tab_elemsize;
 108     int           n0,n1C,n1V,nnn;
 109     real          Y,F,G,H,Fp,Geps,Heps2,epsC,eps2C,epsV,eps2V,VV,FF;
 110     double        isp=0.564189583547756;
 111     real          dvdl_part;
 112
 113     /* fix compiler warnings */
 114     nj1   = 0;
 115     n1C   = n1V   = 0;
 116     epsC  = epsV  = 0;
 117     eps2C = eps2V = 0;
 118
 119     dvdl_coul  = 0;
 120     dvdl_vdw   = 0;
 121
 122     /* Lambda factor for state A, 1-lambda*/
 123     LFC[STATE_A] = 1.0 - lambda_coul;
 124     LFV[STATE_A] = 1.0 - lambda_vdw;
 125
 126     /* Lambda factor for state B, lambda*/
 127     LFC[STATE_B] = lambda_coul;
 128     LFV[STATE_B] = lambda_vdw;
 129
 130     /*derivative of the lambda factor for state A and B */
 131     DLF[STATE_A] = -1;
 132     DLF[STATE_B] = 1;
 133
 134     for (i=0;i<NSTATES;i++)
 135     {
 136         lfac_coul[i]  = (lam_power==2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 137         dlfac_coul[i] = DLF[i]*lam_power/sc_r_power*(lam_power==2 ? (1-LFC[i]) : 1);
 138         lfac_vdw[i]   = (lam_power==2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 139         dlfac_vdw[i]  = DLF[i]*lam_power/sc_r_power*(lam_power==2 ? (1-LFV[i]) : 1);
 140     }
 141     /* precalculate */
 142     sigma2_def = pow(sigma6_def,1.0/3.0);
 143     sigma2_min = pow(sigma6_min,1.0/3.0);
 144
 145     /* Ewald (not PME) table is special (icoul==enbcoulFEWALD) */
 146
 147     do_coultab = (icoul==enbcoulTAB);
 148     do_vdwtab  = (ivdw==enbcoulTAB);
 149
 150     do_tab = do_coultab || do_vdwtab;
 151
 152     /* we always use the combined table here */
 153     tab_elemsize = 12;
 154
 155     for(n=0; (n<nri); n++)
 156     {
 157         is3              = 3*shift[n];
 158         shX              = shiftvec[is3];
 159         shY              = shiftvec[is3+1];
 160         shZ              = shiftvec[is3+2];
 161         nj0              = jindex[n];
 162         nj1              = jindex[n+1];
 163         ii               = iinr[n];
 164         ii3              = 3*ii;
 165         ix               = shX + x[ii3+0];
 166         iy               = shY + x[ii3+1];
 167         iz               = shZ + x[ii3+2];
 168         iqA              = facel*chargeA[ii];
 169         iqB              = facel*chargeB[ii];
 170         ntiA             = 2*ntype*typeA[ii];
 171         ntiB             = 2*ntype*typeB[ii];
 172         vctot            = 0;
 173         vvtot            = 0;
 174         fix              = 0;
 175         fiy              = 0;
 176         fiz              = 0;
 177
 178         for(k=nj0; (k<nj1); k++)
 179         {
 180             jnr              = jjnr[k];
 181             j3               = 3*jnr;
 182             dx               = ix - x[j3];
 183             dy               = iy - x[j3+1];
 184             dz               = iz - x[j3+2];
 185             rsq              = dx*dx+dy*dy+dz*dz;
 186             rinv             = gmx_invsqrt(rsq);
 187             r                = rsq*rinv;
 188             if (sc_r_power == 6.0)
 189             {
 190                 rpm2             = rsq*rsq; /* r4 */
 191                 rp               = rpm2*rsq; /* r6 */
 192             }
 193             else if (sc_r_power == 48.0)
 194             {
 195                 rp               = rsq*rsq*rsq;  /* r6 */
 196                 rp               = rp*rp; /* r12 */
 197                 rp               = rp*rp; /* r24 */
 198                 rp               = rp*rp; /* r48 */
 199                 rpm2             = rp/rsq; /* r46 */
 200             }
 201             else
 202             {
 203                 rp             = pow(r,sc_r_power);  /* not currently supported as input, but can handle it */
 204                 rpm2           = rp/rsq;
 205             }
 206
 207             tj[STATE_A]      = ntiA+2*typeA[jnr];
 208             tj[STATE_B]      = ntiB+2*typeB[jnr];
 209             qq[STATE_A]      = iqA*chargeA[jnr];
 210             qq[STATE_B]      = iqB*chargeB[jnr];
 211
 212             for (i=0;i<NSTATES;i++)
 213             {
 214
 215                 c6[i]              = nbfp[tj[i]];
 216                 c12[i]             = nbfp[tj[i]+1];
 217                 if((c6[i] > 0) && (c12[i] > 0))
 218                 {
 219                     sigma6[i]       = c12[i]/c6[i];
 220                     sigma2[i]       = pow(c12[i]/c6[i],1.0/3.0);
 221                     /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 222                        what data to store externally.  Can't be fixed without larger scale changes, so not 4.6 */
 223                     if (sigma6[i] < sigma6_min) { /* for disappearing coul and vdw with soft core at the same time */
 224                         sigma6[i] = sigma6_min;
 225                         sigma2[i] = sigma2_min;
 226                     }
 227                 }
 228                 else
 229                 {
 230                     sigma6[i]       = sigma6_def;
 231                     sigma2[i]       = sigma2_def;
 232                 }
 233                 if (sc_r_power == 6.0)
 234                 {
 235                     sigma_pow[i]    = sigma6[i];
 236                     sigma_powm2[i]  = sigma6[i]/sigma2[i];
 237                 }
 238                 else if (sc_r_power == 48.0)
 239                 {
 240                     sigma_pow[i]    = sigma6[i]*sigma6[i];   /* sigma^12 */
 241                     sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
 242                     sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
 243                     sigma_powm2[i]  = sigma_pow[i]/sigma2[i];
 244                 }
 245                 else
 246                 {    /* not really supported as input, but in here for testing the general case*/
 247                     sigma_pow[i]    = pow(sigma2[i],sc_r_power/2);
 248                     sigma_powm2[i]  = sigma_pow[i]/(sigma2[i]);
 249                 }
 250             }
 251
 252             /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 253             if ((c12[STATE_A] > 0) && (c12[STATE_B] > 0)) {
 254                 alpha_vdw_eff    = 0;
 255                 alpha_coul_eff   = 0;
 256             } else {
 257                 alpha_vdw_eff    = alpha_vdw;
 258                 alpha_coul_eff   = alpha_coul;
 259             }
 260
 261             for (i=0;i<NSTATES;i++)
 262             {
 263                 FscalC[i]    = 0;
 264                 FscalV[i]    = 0;
 265                 Vcoul[i]     = 0;
 266                 Vvdw[i]      = 0;
 267
 268                 /* Only spend time on A or B state if it is non-zero */
 269                 if( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
 270                 {
 271
 272                     /* this section has to be inside the loop becaue of the dependence on sigma_pow */
 273                     rpinvC         = 1.0/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
 274                     rinvC          = pow(rpinvC,1.0/sc_r_power);
 275                     rC             = 1.0/rinvC;
 276
 277                     rpinvV         = 1.0/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
 278                     rinvV          = pow(rpinvV,1.0/sc_r_power);
 279                     rV             = 1.0/rinvV;
 280
 281                     if (do_tab)
 282                     {
 283                         rtC        = rC*tabscale;
 284                         n0         = rtC;
 285                         epsC       = rtC-n0;
 286                         eps2C      = epsC*epsC;
 287                         n1C        = tab_elemsize*n0;
 288
 289                         rtV        = rV*tabscale;
 290                         n0         = rtV;
 291                         epsV       = rtV-n0;
 292                         eps2V      = epsV*epsV;
 293                         n1V        = tab_elemsize*n0;
 294                     }
 295
 296                     if(icoul==enbcoulOOR || icoul==enbcoulFEWALD)
 297                     {
 298                         /* simple cutoff */
 299                         Vcoul[i]   = qq[i]*rinvC;
 300                         FscalC[i]  = Vcoul[i]*rpinvC;
 301                     }
 302                     else if(icoul==enbcoulRF)
 303                     {
 304                         /* reaction-field */
 305                         krsq       = krf*rC*rC;
 306                         Vcoul[i]   = qq[i]*(rinvC+krsq-crf);
 307                         FscalC[i]  = qq[i]*(rinvC-2.0*krsq)*rpinvC;
 308                     }
 309                     else if (icoul==enbcoulTAB)
 310                     {
 311                         /* non-Ewald tabulated coulomb */
 312                         nnn        = n1C;
 313                         Y          = VFtab[nnn];
 314                         F          = VFtab[nnn+1];
 315                         Geps       = epsC*VFtab[nnn+2];
 316                         Heps2      = eps2C*VFtab[nnn+3];
 317                         Fp         = F+Geps+Heps2;
 318                         VV         = Y+epsC*Fp;
 319                         FF         = Fp+Geps+2.0*Heps2;
 320                         Vcoul[i]   = qq[i]*VV;
 321                         FscalC[i]  = -qq[i]*tabscale*FF*rC*rpinvC;
 322                     }
 323
 324                     if(ivdw==enbvdwLJ)
 325                     {
 326                         /* cutoff LJ */
 327                         if (sc_r_power == 6.0)
 328                         {
 329                             rinv6            = rpinvV;
 330                         }
 331                         else
 332                         {
 333                             rinv6            = pow(rinvV,6.0);
 334                         }
 335                         Vvdw6            = c6[i]*rinv6;
 336                         Vvdw12           = c12[i]*rinv6*rinv6;
 337                         Vvdw[i]          = Vvdw12-Vvdw6;
 338                         FscalV[i]        = (12.0*Vvdw12-6.0*Vvdw6)*rpinvV;
 339                     }
 340                     else if(ivdw==enbvdwTAB)
 341                     {
 342                         /* Table LJ */
 343                         nnn = n1V+4;
 344
 345                         /* dispersion */
 346                         Y          = VFtab[nnn];
 347                         F          = VFtab[nnn+1];
 348                         Geps       = epsV*VFtab[nnn+2];
 349                         Heps2      = eps2V*VFtab[nnn+3];
 350                         Fp         = F+Geps+Heps2;
 351                         VV         = Y+epsV*Fp;
 352                         FF         = Fp+Geps+2.0*Heps2;
 353                         Vvdw[i]   += c6[i]*VV;
 354                         FscalV[i] -= c6[i]*tabscale*FF*rV*rpinvV;
 355
 356                         /* repulsion */
 357                         Y          = VFtab[nnn+4];
 358                         F          = VFtab[nnn+5];
 359                         Geps       = epsV*VFtab[nnn+6];
 360                         Heps2      = eps2V*VFtab[nnn+7];
 361                         Fp         = F+Geps+Heps2;
 362                         VV         = Y+epsV*Fp;
 363                         FF         = Fp+Geps+2.0*Heps2;
 364                         Vvdw[i]   += c12[i]*VV;
 365                         FscalV[i] -= c12[i]*tabscale*FF*rV*rpinvV;
 366                     }
 367                     /* Buckingham vdw free energy not supported for now */
 368                 }
 369             }
 370
 371             Fscal = 0;
 372
 373             if (icoul==enbcoulFEWALD) {
 374                 /* because we compute the softcore normally,
 375                    we have to remove the ewald short range portion. Done outside of
 376                    the states loop because this part doesn't depend on the scaled R */
 377
 378                 if (r != 0)
 379                 {
 380                     VV    = gmx_erf(ewc*r)*rinv;
 381                     FF    = rinv*rinv*(VV - 2.0*ewc*isp*exp(-ewc*ewc*rsq));
 382                 }
 383                 else
 384                 {
 385                     VV    = ewc*2.0/sqrt(M_PI);
 386                     FF    = 0;
 387                 }
 388
 389                 for (i=0;i<NSTATES;i++)
 390                 {
 391                     vctot      -= LFC[i]*qq[i]*VV;
 392                     Fscal      -= LFC[i]*qq[i]*FF;
 393                     dvdl_coul  -= (DLF[i]*qq[i])*VV;
 394                 }
 395             }
 396
 397             /* Assemble A and B states */
 398             for (i=0;i<NSTATES;i++)
 399             {
 400                 vctot         += LFC[i]*Vcoul[i];
 401                 vvtot         += LFV[i]*Vvdw[i];
 402
 403                 Fscal         += LFC[i]*FscalC[i]*rpm2;
 404                 Fscal         += LFV[i]*FscalV[i]*rpm2;
 405
 406                 dvdl_coul     += Vcoul[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*FscalC[i]*sigma_pow[i];
 407                 dvdl_vdw      += Vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*FscalV[i]*sigma_pow[i];
 408             }
 409
 410             if (bDoForces)
 411             {
 412                 tx         = Fscal*dx;
 413                 ty         = Fscal*dy;
 414                 tz         = Fscal*dz;
 415                 fix        = fix + tx;
 416                 fiy        = fiy + ty;
 417                 fiz        = fiz + tz;
 418                 f[j3]      = f[j3]   - tx;
 419                 f[j3+1]    = f[j3+1] - ty;
 420                 f[j3+2]    = f[j3+2] - tz;
 421             }
 422         }
 423
 424         if (bDoForces)
 425         {
 426             f[ii3]         = f[ii3]        + fix;
 427             f[ii3+1]       = f[ii3+1]      + fiy;
 428             f[ii3+2]       = f[ii3+2]      + fiz;
 429             fshift[is3]    = fshift[is3]   + fix;
 430             fshift[is3+1]  = fshift[is3+1] + fiy;
 431             fshift[is3+2]  = fshift[is3+2] + fiz;
 432         }
 433         ggid               = gid[n];
 434         Vc[ggid]           = Vc[ggid] + vctot;
 435         Vv[ggid]           = Vv[ggid] + vvtot;
 436     }
 437
 438     dvdl[efptCOUL]     += dvdl_coul;
 439     dvdl[efptVDW]      += dvdl_vdw;
 440     *outeriter       = nri;
 441     *inneriter       = nj1;
 442 }