?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Redefine the default boolean type to gmx_bool.
[gromacs.git] / src / gmxlib / nonbonded / nb_free_energy.c
blob466c0d4827573e89b0f7e5dd5392d95a02dc92c1
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 }
The ultimate molecular dynamics simulation package