1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
46 gmx_nb_free_energy_kernel(int icoul,
89 int i,j,n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
91 real Fscal,FscalC[NSTATES],FscalV[NSTATES],tx,ty,tz;
92 real Vcoul[NSTATES],Vvdw[NSTATES];
93 real rinv6,r,rt,rtC,rtV;
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;
109 real Y,F,G,H,Fp,Geps,Heps2,epsC,eps2C,epsV,eps2V,VV,FF;
110 double isp=0.564189583547756;
113 /* fix compiler warnings */
122 /* Lambda factor for state A, 1-lambda*/
123 LFC[STATE_A] = 1.0 - lambda_coul;
124 LFV[STATE_A] = 1.0 - lambda_vdw;
126 /* Lambda factor for state B, lambda*/
127 LFC[STATE_B] = lambda_coul;
128 LFV[STATE_B] = lambda_vdw;
130 /*derivative of the lambda factor for state A and B */
134 for (i=0;i<NSTATES;i++)
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);
142 sigma2_def = pow(sigma6_def,1.0/3.0);
143 sigma2_min = pow(sigma6_min,1.0/3.0);
145 /* Ewald (not PME) table is special (icoul==enbcoulFEWALD) */
147 do_coultab = (icoul==enbcoulTAB);
148 do_vdwtab = (ivdw==enbcoulTAB);
150 do_tab = do_coultab || do_vdwtab;
152 /* we always use the combined table here */
155 for(n=0; (n<nri); n++)
159 shY = shiftvec[is3+1];
160 shZ = shiftvec[is3+2];
168 iqA = facel*chargeA[ii];
169 iqB = facel*chargeB[ii];
170 ntiA = 2*ntype*typeA[ii];
171 ntiB = 2*ntype*typeB[ii];
178 for(k=nj0; (k<nj1); k++)
185 rsq = dx*dx+dy*dy+dz*dz;
186 rinv = gmx_invsqrt(rsq);
188 if (sc_r_power == 6.0)
190 rpm2 = rsq*rsq; /* r4 */
191 rp = rpm2*rsq; /* r6 */
193 else if (sc_r_power == 48.0)
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 */
203 rp = pow(r,sc_r_power); /* not currently supported as input, but can handle it */
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];
212 for (i=0;i<NSTATES;i++)
216 c12[i] = nbfp[tj[i]+1];
217 if((c6[i] > 0) && (c12[i] > 0))
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;
230 sigma6[i] = sigma6_def;
231 sigma2[i] = sigma2_def;
233 if (sc_r_power == 6.0)
235 sigma_pow[i] = sigma6[i];
236 sigma_powm2[i] = sigma6[i]/sigma2[i];
238 else if (sc_r_power == 48.0)
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];
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]);
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)) {
257 alpha_vdw_eff = alpha_vdw;
258 alpha_coul_eff = alpha_coul;
261 for (i=0;i<NSTATES;i++)
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) )
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);
277 rpinvV = 1.0/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
278 rinvV = pow(rpinvV,1.0/sc_r_power);
287 n1C = tab_elemsize*n0;
293 n1V = tab_elemsize*n0;
296 if(icoul==enbcoulOOR || icoul==enbcoulFEWALD)
299 Vcoul[i] = qq[i]*rinvC;
300 FscalC[i] = Vcoul[i]*rpinvC;
302 else if(icoul==enbcoulRF)
306 Vcoul[i] = qq[i]*(rinvC+krsq-crf);
307 FscalC[i] = qq[i]*(rinvC-2.0*krsq)*rpinvC;
309 else if (icoul==enbcoulTAB)
311 /* non-Ewald tabulated coulomb */
315 Geps = epsC*VFtab[nnn+2];
316 Heps2 = eps2C*VFtab[nnn+3];
319 FF = Fp+Geps+2.0*Heps2;
321 FscalC[i] = -qq[i]*tabscale*FF*rC*rpinvC;
327 if (sc_r_power == 6.0)
333 rinv6 = pow(rinvV,6.0);
336 Vvdw12 = c12[i]*rinv6*rinv6;
337 Vvdw[i] = Vvdw12-Vvdw6;
338 FscalV[i] = (12.0*Vvdw12-6.0*Vvdw6)*rpinvV;
340 else if(ivdw==enbvdwTAB)
348 Geps = epsV*VFtab[nnn+2];
349 Heps2 = eps2V*VFtab[nnn+3];
352 FF = Fp+Geps+2.0*Heps2;
354 FscalV[i] -= c6[i]*tabscale*FF*rV*rpinvV;
359 Geps = epsV*VFtab[nnn+6];
360 Heps2 = eps2V*VFtab[nnn+7];
363 FF = Fp+Geps+2.0*Heps2;
364 Vvdw[i] += c12[i]*VV;
365 FscalV[i] -= c12[i]*tabscale*FF*rV*rpinvV;
367 /* Buckingham vdw free energy not supported for now */
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 */
380 VV = gmx_erf(ewc*r)*rinv;
381 FF = rinv*rinv*(VV - 2.0*ewc*isp*exp(-ewc*ewc*rsq));
385 VV = ewc*2.0/sqrt(M_PI);
389 for (i=0;i<NSTATES;i++)
391 vctot -= LFC[i]*qq[i]*VV;
392 Fscal -= LFC[i]*qq[i]*FF;
393 dvdl_coul -= (DLF[i]*qq[i])*VV;
397 /* Assemble A and B states */
398 for (i=0;i<NSTATES;i++)
400 vctot += LFC[i]*Vcoul[i];
401 vvtot += LFV[i]*Vvdw[i];
403 Fscal += LFC[i]*FscalC[i]*rpm2;
404 Fscal += LFV[i]*FscalV[i]*rpm2;
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];
419 f[j3+1] = f[j3+1] - ty;
420 f[j3+2] = f[j3+2] - tz;
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;
434 Vc[ggid] = Vc[ggid] + vctot;
435 Vv[ggid] = Vv[ggid] + vvtot;
438 dvdl[efptCOUL] += dvdl_coul;
439 dvdl[efptVDW] += dvdl_vdw;