File: | gromacs/gmxlib/nonbonded/nb_free_energy.c |
Location: | line 247, column 5 |
Description: | Value stored to 'nj1' is never read |
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 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <math.h> |
42 | |
43 | #include "gromacs/math/vec.h" |
44 | #include "typedefs.h" |
45 | #include "nonbonded.h" |
46 | #include "nb_kernel.h" |
47 | #include "nrnb.h" |
48 | #include "macros.h" |
49 | #include "nb_free_energy.h" |
50 | |
51 | #include "gromacs/utility/fatalerror.h" |
52 | |
53 | void |
54 | gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict__restrict nlist, |
55 | rvec * gmx_restrict__restrict xx, |
56 | rvec * gmx_restrict__restrict ff, |
57 | t_forcerec * gmx_restrict__restrict fr, |
58 | const t_mdatoms * gmx_restrict__restrict mdatoms, |
59 | nb_kernel_data_t * gmx_restrict__restrict kernel_data, |
60 | t_nrnb * gmx_restrict__restrict nrnb) |
61 | { |
62 | |
63 | #define STATE_A0 0 |
64 | #define STATE_B1 1 |
65 | #define NSTATES2 2 |
66 | int i, j, n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid; |
67 | real shX, shY, shZ; |
68 | real Fscal, FscalC[NSTATES2], FscalV[NSTATES2], tx, ty, tz; |
69 | real Vcoul[NSTATES2], Vvdw[NSTATES2]; |
70 | real rinv6, r, rt, rtC, rtV; |
71 | real iqA, iqB; |
72 | real qq[NSTATES2], vctot, krsq; |
73 | int ntiA, ntiB, tj[NSTATES2]; |
74 | real Vvdw6, Vvdw12, vvtot; |
75 | real ix, iy, iz, fix, fiy, fiz; |
76 | real dx, dy, dz, rsq, rinv; |
77 | real c6[NSTATES2], c12[NSTATES2], c6grid[NSTATES2]; |
78 | real LFC[NSTATES2], LFV[NSTATES2], DLF[NSTATES2]; |
79 | double dvdl_coul, dvdl_vdw; |
80 | real lfac_coul[NSTATES2], dlfac_coul[NSTATES2], lfac_vdw[NSTATES2], dlfac_vdw[NSTATES2]; |
81 | real sigma6[NSTATES2], alpha_vdw_eff, alpha_coul_eff, sigma2_def, sigma2_min; |
82 | real rp, rpm2, rC, rV, rinvC, rpinvC, rinvV, rpinvV; |
83 | real sigma2[NSTATES2], sigma_pow[NSTATES2], sigma_powm2[NSTATES2], rs, rs2; |
84 | int do_tab, tab_elemsize; |
85 | int n0, n1C, n1V, nnn; |
86 | real Y, F, G, H, Fp, Geps, Heps2, epsC, eps2C, epsV, eps2V, VV, FF; |
87 | int icoul, ivdw; |
88 | int nri; |
89 | const int * iinr; |
90 | const int * jindex; |
91 | const int * jjnr; |
92 | const int * shift; |
93 | const int * gid; |
94 | const int * typeA; |
95 | const int * typeB; |
96 | int ntype; |
97 | const real * shiftvec; |
98 | real dvdl_part; |
99 | real * fshift; |
100 | real tabscale = 0; |
101 | const real * VFtab = NULL((void*)0); |
102 | const real * x; |
103 | real * f; |
104 | real facel, krf, crf; |
105 | const real * chargeA; |
106 | const real * chargeB; |
107 | real sigma6_min, sigma6_def, lam_power, sc_power, sc_r_power; |
108 | real alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc_lj; |
109 | const real * nbfp, *nbfp_grid; |
110 | real * dvdl; |
111 | real * Vv; |
112 | real * Vc; |
113 | gmx_bool bDoForces, bDoShiftForces, bDoPotential; |
114 | real rcoulomb, sh_ewald; |
115 | real rvdw, sh_invrc6; |
116 | gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoffAll, bEwald; |
117 | real rcutoff_max2; |
118 | real rcutoff, rcutoff2, rswitch, d, d2, swV3, swV4, swV5, swF2, swF3, swF4, sw, dsw, rinvcorr; |
119 | const real * tab_ewald_F; |
120 | const real * tab_ewald_V; |
121 | const real * tab_ewald_F_lj; |
122 | const real * tab_ewald_V_lj; |
123 | real tab_ewald_scale, tab_ewald_halfsp; |
124 | |
125 | x = xx[0]; |
126 | f = ff[0]; |
127 | |
128 | fshift = fr->fshift[0]; |
129 | |
130 | nri = nlist->nri; |
131 | iinr = nlist->iinr; |
132 | jindex = nlist->jindex; |
133 | jjnr = nlist->jjnr; |
134 | icoul = nlist->ielec; |
135 | ivdw = nlist->ivdw; |
136 | shift = nlist->shift; |
137 | gid = nlist->gid; |
138 | |
139 | shiftvec = fr->shift_vec[0]; |
140 | chargeA = mdatoms->chargeA; |
141 | chargeB = mdatoms->chargeB; |
142 | facel = fr->epsfac; |
143 | krf = fr->k_rf; |
144 | crf = fr->c_rf; |
145 | ewc_lj = fr->ewaldcoeff_lj; |
146 | Vc = kernel_data->energygrp_elec; |
147 | typeA = mdatoms->typeA; |
148 | typeB = mdatoms->typeB; |
149 | ntype = fr->ntype; |
150 | nbfp = fr->nbfp; |
151 | nbfp_grid = fr->ljpme_c6grid; |
152 | Vv = kernel_data->energygrp_vdw; |
153 | lambda_coul = kernel_data->lambda[efptCOUL]; |
154 | lambda_vdw = kernel_data->lambda[efptVDW]; |
155 | dvdl = kernel_data->dvdl; |
156 | alpha_coul = fr->sc_alphacoul; |
157 | alpha_vdw = fr->sc_alphavdw; |
158 | lam_power = fr->sc_power; |
159 | sc_r_power = fr->sc_r_power; |
160 | sigma6_def = fr->sc_sigma6_def; |
161 | sigma6_min = fr->sc_sigma6_min; |
162 | bDoForces = kernel_data->flags & GMX_NONBONDED_DO_FORCE(1<<1); |
163 | bDoShiftForces = kernel_data->flags & GMX_NONBONDED_DO_SHIFTFORCE(1<<2); |
164 | bDoPotential = kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL(1<<4); |
165 | |
166 | rcoulomb = fr->rcoulomb; |
167 | sh_ewald = fr->ic->sh_ewald; |
168 | rvdw = fr->rvdw; |
169 | sh_invrc6 = fr->ic->sh_invrc6; |
170 | |
171 | /* Ewald (PME) reciprocal force and energy quadratic spline tables */ |
172 | tab_ewald_F = fr->ic->tabq_coul_F; |
173 | tab_ewald_V = fr->ic->tabq_coul_V; |
174 | tab_ewald_scale = fr->ic->tabq_scale; |
175 | tab_ewald_F_lj = fr->ic->tabq_vdw_F; |
176 | tab_ewald_V_lj = fr->ic->tabq_vdw_V; |
177 | tab_ewald_halfsp = 0.5/tab_ewald_scale; |
178 | |
179 | if (fr->coulomb_modifier == eintmodPOTSWITCH || fr->vdw_modifier == eintmodPOTSWITCH) |
180 | { |
181 | rcutoff = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb : fr->rvdw; |
182 | rcutoff2 = rcutoff*rcutoff; |
183 | rswitch = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb_switch : fr->rvdw_switch; |
184 | d = rcutoff-rswitch; |
185 | swV3 = -10.0/(d*d*d); |
186 | swV4 = 15.0/(d*d*d*d); |
187 | swV5 = -6.0/(d*d*d*d*d); |
188 | swF2 = -30.0/(d*d*d); |
189 | swF3 = 60.0/(d*d*d*d); |
190 | swF4 = -30.0/(d*d*d*d*d); |
191 | } |
192 | else |
193 | { |
194 | /* Stupid compilers dont realize these variables will not be used */ |
195 | rswitch = 0.0; |
196 | swV3 = 0.0; |
197 | swV4 = 0.0; |
198 | swV5 = 0.0; |
199 | swF2 = 0.0; |
200 | swF3 = 0.0; |
201 | swF4 = 0.0; |
202 | } |
203 | |
204 | if (fr->cutoff_scheme == ecutsVERLET) |
205 | { |
206 | const interaction_const_t *ic; |
207 | |
208 | ic = fr->ic; |
209 | if (EVDW_PME(ic->vdwtype)((ic->vdwtype) == evdwPME)) |
210 | { |
211 | ivdw = GMX_NBKERNEL_VDW_LJEWALD; |
212 | } |
213 | else |
214 | { |
215 | ivdw = GMX_NBKERNEL_VDW_LENNARDJONES; |
216 | } |
217 | |
218 | if (ic->eeltype == eelCUT || EEL_RF(ic->eeltype)((ic->eeltype) == eelRF || (ic->eeltype) == eelGRF || ( ic->eeltype) == eelRF_NEC || (ic->eeltype) == eelRF_ZERO )) |
219 | { |
220 | icoul = GMX_NBKERNEL_ELEC_REACTIONFIELD; |
221 | } |
222 | else if (EEL_PME_EWALD(ic->eeltype)(((ic->eeltype) == eelPME || (ic->eeltype) == eelPMESWITCH || (ic->eeltype) == eelPMEUSER || (ic->eeltype) == eelPMEUSERSWITCH || (ic->eeltype) == eelP3M_AD) || (ic->eeltype) == eelEWALD )) |
223 | { |
224 | icoul = GMX_NBKERNEL_ELEC_EWALD; |
225 | } |
226 | else |
227 | { |
228 | gmx_incons("Unsupported eeltype with Verlet and free-energy")_gmx_error("incons", "Unsupported eeltype with Verlet and free-energy" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/nonbonded/nb_free_energy.c" , 228); |
229 | } |
230 | |
231 | bExactElecCutoff = TRUE1; |
232 | bExactVdwCutoff = TRUE1; |
233 | } |
234 | else |
235 | { |
236 | bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO; |
237 | bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE); |
238 | } |
239 | |
240 | bExactCutoffAll = (bExactElecCutoff && bExactVdwCutoff); |
241 | rcutoff_max2 = max(fr->rcoulomb, fr->rvdw)(((fr->rcoulomb) > (fr->rvdw)) ? (fr->rcoulomb) : (fr->rvdw) ); |
242 | rcutoff_max2 = rcutoff_max2*rcutoff_max2; |
243 | |
244 | bEwald = (icoul == GMX_NBKERNEL_ELEC_EWALD); |
245 | |
246 | /* fix compiler warnings */ |
247 | nj1 = 0; |
Value stored to 'nj1' is never read | |
248 | n1C = n1V = 0; |
249 | epsC = epsV = 0; |
250 | eps2C = eps2V = 0; |
251 | |
252 | dvdl_coul = 0; |
253 | dvdl_vdw = 0; |
254 | |
255 | /* Lambda factor for state A, 1-lambda*/ |
256 | LFC[STATE_A0] = 1.0 - lambda_coul; |
257 | LFV[STATE_A0] = 1.0 - lambda_vdw; |
258 | |
259 | /* Lambda factor for state B, lambda*/ |
260 | LFC[STATE_B1] = lambda_coul; |
261 | LFV[STATE_B1] = lambda_vdw; |
262 | |
263 | /*derivative of the lambda factor for state A and B */ |
264 | DLF[STATE_A0] = -1; |
265 | DLF[STATE_B1] = 1; |
266 | |
267 | for (i = 0; i < NSTATES2; i++) |
268 | { |
269 | lfac_coul[i] = (lam_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i])); |
270 | dlfac_coul[i] = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFC[i]) : 1); |
271 | lfac_vdw[i] = (lam_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i])); |
272 | dlfac_vdw[i] = DLF[i]*lam_power/sc_r_power*(lam_power == 2 ? (1-LFV[i]) : 1); |
273 | } |
274 | /* precalculate */ |
275 | sigma2_def = pow(sigma6_def, 1.0/3.0); |
276 | sigma2_min = pow(sigma6_min, 1.0/3.0); |
277 | |
278 | /* Ewald (not PME) table is special (icoul==enbcoulFEWALD) */ |
279 | |
280 | do_tab = (icoul == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE || |
281 | ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE); |
282 | if (do_tab) |
283 | { |
284 | tabscale = kernel_data->table_elec_vdw->scale; |
285 | VFtab = kernel_data->table_elec_vdw->data; |
286 | /* we always use the combined table here */ |
287 | tab_elemsize = 12; |
288 | } |
289 | |
290 | for (n = 0; (n < nri); n++) |
291 | { |
292 | int npair_within_cutoff; |
293 | |
294 | npair_within_cutoff = 0; |
295 | |
296 | is3 = 3*shift[n]; |
297 | shX = shiftvec[is3]; |
298 | shY = shiftvec[is3+1]; |
299 | shZ = shiftvec[is3+2]; |
300 | nj0 = jindex[n]; |
301 | nj1 = jindex[n+1]; |
302 | ii = iinr[n]; |
303 | ii3 = 3*ii; |
304 | ix = shX + x[ii3+0]; |
305 | iy = shY + x[ii3+1]; |
306 | iz = shZ + x[ii3+2]; |
307 | iqA = facel*chargeA[ii]; |
308 | iqB = facel*chargeB[ii]; |
309 | ntiA = 2*ntype*typeA[ii]; |
310 | ntiB = 2*ntype*typeB[ii]; |
311 | vctot = 0; |
312 | vvtot = 0; |
313 | fix = 0; |
314 | fiy = 0; |
315 | fiz = 0; |
316 | |
317 | for (k = nj0; (k < nj1); k++) |
318 | { |
319 | jnr = jjnr[k]; |
320 | j3 = 3*jnr; |
321 | dx = ix - x[j3]; |
322 | dy = iy - x[j3+1]; |
323 | dz = iz - x[j3+2]; |
324 | rsq = dx*dx + dy*dy + dz*dz; |
325 | |
326 | if (bExactCutoffAll && rsq >= rcutoff_max2) |
327 | { |
328 | /* We save significant time by skipping all code below. |
329 | * Note that with soft-core interactions, the actual cut-off |
330 | * check might be different. But since the soft-core distance |
331 | * is always larger than r, checking on r here is safe. |
332 | */ |
333 | continue; |
334 | } |
335 | npair_within_cutoff++; |
336 | |
337 | if (rsq > 0) |
338 | { |
339 | rinv = gmx_invsqrt(rsq)gmx_software_invsqrt(rsq); |
340 | r = rsq*rinv; |
341 | } |
342 | else |
343 | { |
344 | /* The force at r=0 is zero, because of symmetry. |
345 | * But note that the potential is in general non-zero, |
346 | * since the soft-cored r will be non-zero. |
347 | */ |
348 | rinv = 0; |
349 | r = 0; |
350 | } |
351 | |
352 | if (sc_r_power == 6.0) |
353 | { |
354 | rpm2 = rsq*rsq; /* r4 */ |
355 | rp = rpm2*rsq; /* r6 */ |
356 | } |
357 | else if (sc_r_power == 48.0) |
358 | { |
359 | rp = rsq*rsq*rsq; /* r6 */ |
360 | rp = rp*rp; /* r12 */ |
361 | rp = rp*rp; /* r24 */ |
362 | rp = rp*rp; /* r48 */ |
363 | rpm2 = rp/rsq; /* r46 */ |
364 | } |
365 | else |
366 | { |
367 | rp = pow(r, sc_r_power); /* not currently supported as input, but can handle it */ |
368 | rpm2 = rp/rsq; |
369 | } |
370 | |
371 | Fscal = 0; |
372 | |
373 | qq[STATE_A0] = iqA*chargeA[jnr]; |
374 | qq[STATE_B1] = iqB*chargeB[jnr]; |
375 | |
376 | tj[STATE_A0] = ntiA+2*typeA[jnr]; |
377 | tj[STATE_B1] = ntiB+2*typeB[jnr]; |
378 | |
379 | if (ivdw == GMX_NBKERNEL_VDW_LJEWALD) |
380 | { |
381 | c6grid[STATE_A0] = nbfp_grid[tj[STATE_A0]]; |
382 | c6grid[STATE_B1] = nbfp_grid[tj[STATE_B1]]; |
383 | } |
384 | |
385 | if (nlist->excl_fep == NULL((void*)0) || nlist->excl_fep[k]) |
386 | { |
387 | c6[STATE_A0] = nbfp[tj[STATE_A0]]; |
388 | c6[STATE_B1] = nbfp[tj[STATE_B1]]; |
389 | |
390 | for (i = 0; i < NSTATES2; i++) |
391 | { |
392 | c12[i] = nbfp[tj[i]+1]; |
393 | if ((c6[i] > 0) && (c12[i] > 0)) |
394 | { |
395 | /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */ |
396 | sigma6[i] = 0.5*c12[i]/c6[i]; |
397 | sigma2[i] = pow(sigma6[i], 1.0/3.0); |
398 | /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on |
399 | what data to store externally. Can't be fixed without larger scale changes, so not 4.6 */ |
400 | if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */ |
401 | { |
402 | sigma6[i] = sigma6_min; |
403 | sigma2[i] = sigma2_min; |
404 | } |
405 | } |
406 | else |
407 | { |
408 | sigma6[i] = sigma6_def; |
409 | sigma2[i] = sigma2_def; |
410 | } |
411 | if (sc_r_power == 6.0) |
412 | { |
413 | sigma_pow[i] = sigma6[i]; |
414 | sigma_powm2[i] = sigma6[i]/sigma2[i]; |
415 | } |
416 | else if (sc_r_power == 48.0) |
417 | { |
418 | sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */ |
419 | sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */ |
420 | sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */ |
421 | sigma_powm2[i] = sigma_pow[i]/sigma2[i]; |
422 | } |
423 | else |
424 | { /* not really supported as input, but in here for testing the general case*/ |
425 | sigma_pow[i] = pow(sigma2[i], sc_r_power/2); |
426 | sigma_powm2[i] = sigma_pow[i]/(sigma2[i]); |
427 | } |
428 | } |
429 | |
430 | /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ |
431 | if ((c12[STATE_A0] > 0) && (c12[STATE_B1] > 0)) |
432 | { |
433 | alpha_vdw_eff = 0; |
434 | alpha_coul_eff = 0; |
435 | } |
436 | else |
437 | { |
438 | alpha_vdw_eff = alpha_vdw; |
439 | alpha_coul_eff = alpha_coul; |
440 | } |
441 | |
442 | for (i = 0; i < NSTATES2; i++) |
443 | { |
444 | FscalC[i] = 0; |
445 | FscalV[i] = 0; |
446 | Vcoul[i] = 0; |
447 | Vvdw[i] = 0; |
448 | |
449 | /* Only spend time on A or B state if it is non-zero */ |
450 | if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) ) |
451 | { |
452 | /* this section has to be inside the loop because of the dependence on sigma_pow */ |
453 | rpinvC = 1.0/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp); |
454 | rinvC = pow(rpinvC, 1.0/sc_r_power); |
455 | rC = 1.0/rinvC; |
456 | |
457 | rpinvV = 1.0/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp); |
458 | rinvV = pow(rpinvV, 1.0/sc_r_power); |
459 | rV = 1.0/rinvV; |
460 | |
461 | if (do_tab) |
462 | { |
463 | rtC = rC*tabscale; |
464 | n0 = rtC; |
465 | epsC = rtC-n0; |
466 | eps2C = epsC*epsC; |
467 | n1C = tab_elemsize*n0; |
468 | |
469 | rtV = rV*tabscale; |
470 | n0 = rtV; |
471 | epsV = rtV-n0; |
472 | eps2V = epsV*epsV; |
473 | n1V = tab_elemsize*n0; |
474 | } |
475 | |
476 | /* With Ewald and soft-core we should put the cut-off on r, |
477 | * not on the soft-cored rC, as the real-space and |
478 | * reciprocal space contributions should (almost) cancel. |
479 | */ |
480 | if (qq[i] != 0 && |
481 | !(bExactElecCutoff && |
482 | ((!bEwald && rC >= rcoulomb) || |
483 | (bEwald && r >= rcoulomb)))) |
484 | { |
485 | switch (icoul) |
486 | { |
487 | case GMX_NBKERNEL_ELEC_COULOMB: |
488 | /* simple cutoff */ |
489 | Vcoul[i] = qq[i]*rinvC; |
490 | FscalC[i] = Vcoul[i]; |
491 | break; |
492 | |
493 | case GMX_NBKERNEL_ELEC_EWALD: |
494 | /* Ewald FEP is done only on the 1/r part */ |
495 | Vcoul[i] = qq[i]*(rinvC - sh_ewald); |
496 | FscalC[i] = Vcoul[i]; |
497 | break; |
498 | |
499 | case GMX_NBKERNEL_ELEC_REACTIONFIELD: |
500 | /* reaction-field */ |
501 | Vcoul[i] = qq[i]*(rinvC + krf*rC*rC-crf); |
502 | FscalC[i] = qq[i]*(rinvC - 2.0*krf*rC*rC); |
503 | break; |
504 | |
505 | case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE: |
506 | /* non-Ewald tabulated coulomb */ |
507 | nnn = n1C; |
508 | Y = VFtab[nnn]; |
509 | F = VFtab[nnn+1]; |
510 | Geps = epsC*VFtab[nnn+2]; |
511 | Heps2 = eps2C*VFtab[nnn+3]; |
512 | Fp = F+Geps+Heps2; |
513 | VV = Y+epsC*Fp; |
514 | FF = Fp+Geps+2.0*Heps2; |
515 | Vcoul[i] = qq[i]*VV; |
516 | FscalC[i] = -qq[i]*tabscale*FF*rC; |
517 | break; |
518 | |
519 | case GMX_NBKERNEL_ELEC_GENERALIZEDBORN: |
520 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/nonbonded/nb_free_energy.c" , 520, "Free energy and GB not implemented.\n"); |
521 | break; |
522 | |
523 | case GMX_NBKERNEL_ELEC_NONE: |
524 | FscalC[i] = 0.0; |
525 | Vcoul[i] = 0.0; |
526 | break; |
527 | |
528 | default: |
529 | gmx_incons("Invalid icoul in free energy kernel")_gmx_error("incons", "Invalid icoul in free energy kernel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/nonbonded/nb_free_energy.c" , 529); |
530 | break; |
531 | } |
532 | |
533 | if (fr->coulomb_modifier == eintmodPOTSWITCH) |
534 | { |
535 | d = rC-rswitch; |
536 | d = (d > 0.0) ? d : 0.0; |
537 | d2 = d*d; |
538 | sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5)); |
539 | dsw = d2*(swF2+d*(swF3+d*swF4)); |
540 | |
541 | Vcoul[i] *= sw; |
542 | FscalC[i] = FscalC[i]*sw + Vcoul[i]*dsw; |
543 | } |
544 | } |
545 | |
546 | if ((c6[i] != 0 || c12[i] != 0) && |
547 | !(bExactVdwCutoff && |
548 | ((ivdw != GMX_NBKERNEL_VDW_LJEWALD && rV >= rvdw) || |
549 | (ivdw == GMX_NBKERNEL_VDW_LJEWALD && r >= rvdw)))) |
550 | { |
551 | switch (ivdw) |
552 | { |
553 | case GMX_NBKERNEL_VDW_LENNARDJONES: |
554 | case GMX_NBKERNEL_VDW_LJEWALD: |
555 | /* cutoff LJ */ |
556 | if (sc_r_power == 6.0) |
557 | { |
558 | rinv6 = rpinvV; |
559 | } |
560 | else |
561 | { |
562 | rinv6 = pow(rinvV, 6.0); |
563 | } |
564 | Vvdw6 = c6[i]*rinv6; |
565 | Vvdw12 = c12[i]*rinv6*rinv6; |
566 | if (fr->vdw_modifier == eintmodPOTSHIFT) |
567 | { |
568 | Vvdw[i] = ( (Vvdw12-c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0) |
569 | -(Vvdw6-c6[i]*sh_invrc6)*(1.0/6.0)); |
570 | } |
571 | else |
572 | { |
573 | Vvdw[i] = Vvdw12*(1.0/12.0) - Vvdw6*(1.0/6.0); |
574 | } |
575 | FscalV[i] = Vvdw12 - Vvdw6; |
576 | break; |
577 | |
578 | case GMX_NBKERNEL_VDW_BUCKINGHAM: |
579 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/nonbonded/nb_free_energy.c" , 579, "Buckingham free energy not supported."); |
580 | break; |
581 | |
582 | case GMX_NBKERNEL_VDW_CUBICSPLINETABLE: |
583 | /* Table LJ */ |
584 | nnn = n1V+4; |
585 | /* dispersion */ |
586 | Y = VFtab[nnn]; |
587 | F = VFtab[nnn+1]; |
588 | Geps = epsV*VFtab[nnn+2]; |
589 | Heps2 = eps2V*VFtab[nnn+3]; |
590 | Fp = F+Geps+Heps2; |
591 | VV = Y+epsV*Fp; |
592 | FF = Fp+Geps+2.0*Heps2; |
593 | Vvdw[i] += c6[i]*VV; |
594 | FscalV[i] -= c6[i]*tabscale*FF*rV; |
595 | |
596 | /* repulsion */ |
597 | Y = VFtab[nnn+4]; |
598 | F = VFtab[nnn+5]; |
599 | Geps = epsV*VFtab[nnn+6]; |
600 | Heps2 = eps2V*VFtab[nnn+7]; |
601 | Fp = F+Geps+Heps2; |
602 | VV = Y+epsV*Fp; |
603 | FF = Fp+Geps+2.0*Heps2; |
604 | Vvdw[i] += c12[i]*VV; |
605 | FscalV[i] -= c12[i]*tabscale*FF*rV; |
606 | break; |
607 | |
608 | case GMX_NBKERNEL_VDW_NONE: |
609 | Vvdw[i] = 0.0; |
610 | FscalV[i] = 0.0; |
611 | break; |
612 | |
613 | default: |
614 | gmx_incons("Invalid ivdw in free energy kernel")_gmx_error("incons", "Invalid ivdw in free energy kernel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/nonbonded/nb_free_energy.c" , 614); |
615 | break; |
616 | } |
617 | |
618 | if (fr->vdw_modifier == eintmodPOTSWITCH) |
619 | { |
620 | d = rV-rswitch; |
621 | d = (d > 0.0) ? d : 0.0; |
622 | d2 = d*d; |
623 | sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5)); |
624 | dsw = d2*(swF2+d*(swF3+d*swF4)); |
625 | |
626 | Vvdw[i] *= sw; |
627 | FscalV[i] = FscalV[i]*sw + Vvdw[i]*dsw; |
628 | |
629 | FscalV[i] = (rV < rvdw) ? FscalV[i] : 0.0; |
630 | Vvdw[i] = (rV < rvdw) ? Vvdw[i] : 0.0; |
631 | } |
632 | } |
633 | |
634 | /* FscalC (and FscalV) now contain: dV/drC * rC |
635 | * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p |
636 | * Further down we first multiply by r^p-2 and then by |
637 | * the vector r, which in total gives: dV/drC * (r/rC)^1-p |
638 | */ |
639 | FscalC[i] *= rpinvC; |
640 | FscalV[i] *= rpinvV; |
641 | } |
642 | } |
643 | |
644 | /* Assemble A and B states */ |
645 | for (i = 0; i < NSTATES2; i++) |
646 | { |
647 | vctot += LFC[i]*Vcoul[i]; |
648 | vvtot += LFV[i]*Vvdw[i]; |
649 | |
650 | Fscal += LFC[i]*FscalC[i]*rpm2; |
651 | Fscal += LFV[i]*FscalV[i]*rpm2; |
652 | |
653 | dvdl_coul += Vcoul[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*FscalC[i]*sigma_pow[i]; |
654 | dvdl_vdw += Vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*FscalV[i]*sigma_pow[i]; |
655 | } |
656 | } |
657 | else if (icoul == GMX_NBKERNEL_ELEC_REACTIONFIELD) |
658 | { |
659 | /* For excluded pairs, which are only in this pair list when |
660 | * using the Verlet scheme, we don't use soft-core. |
661 | * The group scheme also doesn't soft-core for these. |
662 | * As there is no singularity, there is no need for soft-core. |
663 | */ |
664 | VV = krf*rsq - crf; |
665 | FF = -2.0*krf; |
666 | |
667 | if (ii == jnr) |
668 | { |
669 | VV *= 0.5; |
670 | } |
671 | |
672 | for (i = 0; i < NSTATES2; i++) |
673 | { |
674 | vctot += LFC[i]*qq[i]*VV; |
675 | Fscal += LFC[i]*qq[i]*FF; |
676 | dvdl_coul += DLF[i]*qq[i]*VV; |
677 | } |
678 | } |
679 | |
680 | if (icoul == GMX_NBKERNEL_ELEC_EWALD && |
681 | !(bExactElecCutoff && r >= rcoulomb)) |
682 | { |
683 | /* Because we compute the soft-core normally, |
684 | * we have to remove the Ewald short range portion. |
685 | * Done outside of the states loop because this part |
686 | * doesn't depend on the scaled R. |
687 | */ |
688 | real rs, frac, f_lr; |
689 | int ri; |
690 | |
691 | rs = rsq*rinv*tab_ewald_scale; |
692 | ri = (int)rs; |
693 | frac = rs - ri; |
694 | f_lr = (1 - frac)*tab_ewald_F[ri] + frac*tab_ewald_F[ri+1]; |
695 | FF = f_lr*rinv; |
696 | VV = tab_ewald_V[ri] - tab_ewald_halfsp*frac*(tab_ewald_F[ri] + f_lr); |
697 | |
698 | if (ii == jnr) |
699 | { |
700 | VV *= 0.5; |
701 | } |
702 | |
703 | for (i = 0; i < NSTATES2; i++) |
704 | { |
705 | vctot -= LFC[i]*qq[i]*VV; |
706 | Fscal -= LFC[i]*qq[i]*FF; |
707 | dvdl_coul -= (DLF[i]*qq[i])*VV; |
708 | } |
709 | } |
710 | |
711 | if (ivdw == GMX_NBKERNEL_VDW_LJEWALD && |
712 | !(bExactVdwCutoff && r >= rvdw)) |
713 | { |
714 | real rs, frac, f_lr; |
715 | int ri; |
716 | |
717 | rs = rsq*rinv*tab_ewald_scale; |
718 | ri = (int)rs; |
719 | frac = rs - ri; |
720 | f_lr = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1]; |
721 | FF = f_lr*rinv; |
722 | VV = tab_ewald_V_lj[ri] - tab_ewald_halfsp*frac*(tab_ewald_F_lj[ri] + f_lr); |
723 | for (i = 0; i < NSTATES2; i++) |
724 | { |
725 | vvtot += LFV[i]*c6grid[i]*VV*(1.0/6.0); |
726 | Fscal += LFV[i]*c6grid[i]*FF*(1.0/6.0); |
727 | dvdl_vdw += (DLF[i]*c6grid[i])*VV*(1.0/6.0); |
728 | } |
729 | |
730 | } |
731 | |
732 | if (bDoForces) |
733 | { |
734 | tx = Fscal*dx; |
735 | ty = Fscal*dy; |
736 | tz = Fscal*dz; |
737 | fix = fix + tx; |
738 | fiy = fiy + ty; |
739 | fiz = fiz + tz; |
740 | /* OpenMP atomics are expensive, but this kernels is also |
741 | * expensive, so we can take this hit, instead of using |
742 | * thread-local output buffers and extra reduction. |
743 | */ |
744 | #pragma omp atomic |
745 | f[j3] -= tx; |
746 | #pragma omp atomic |
747 | f[j3+1] -= ty; |
748 | #pragma omp atomic |
749 | f[j3+2] -= tz; |
750 | } |
751 | } |
752 | |
753 | /* The atomics below are expensive with many OpenMP threads. |
754 | * Here unperturbed i-particles will usually only have a few |
755 | * (perturbed) j-particles in the list. Thus with a buffered list |
756 | * we can skip a significant number of i-reductions with a check. |
757 | */ |
758 | if (npair_within_cutoff > 0) |
759 | { |
760 | if (bDoForces) |
761 | { |
762 | #pragma omp atomic |
763 | f[ii3] += fix; |
764 | #pragma omp atomic |
765 | f[ii3+1] += fiy; |
766 | #pragma omp atomic |
767 | f[ii3+2] += fiz; |
768 | } |
769 | if (bDoShiftForces) |
770 | { |
771 | #pragma omp atomic |
772 | fshift[is3] += fix; |
773 | #pragma omp atomic |
774 | fshift[is3+1] += fiy; |
775 | #pragma omp atomic |
776 | fshift[is3+2] += fiz; |
777 | } |
778 | if (bDoPotential) |
779 | { |
780 | ggid = gid[n]; |
781 | #pragma omp atomic |
782 | Vc[ggid] += vctot; |
783 | #pragma omp atomic |
784 | Vv[ggid] += vvtot; |
785 | } |
786 | } |
787 | } |
788 | |
789 | #pragma omp atomic |
790 | dvdl[efptCOUL] += dvdl_coul; |
791 | #pragma omp atomic |
792 | dvdl[efptVDW] += dvdl_vdw; |
793 | |
794 | /* Estimate flops, average for free energy stuff: |
795 | * 12 flops per outer iteration |
796 | * 150 flops per inner iteration |
797 | */ |
798 | inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150)(nrnb)->n[eNR_NBKERNEL_FREE_ENERGY] += nlist->nri*12 + nlist ->jindex[n]*150; |
799 | } |
800 | |
801 | real |
802 | nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, real alpha_vdw, |
803 | real tabscale, real *vftab, |
804 | real qqA, real c6A, real c12A, real qqB, real c6B, real c12B, |
805 | real LFC[2], real LFV[2], real DLF[2], |
806 | real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], real dlfac_vdw[2], |
807 | real sigma6_def, real sigma6_min, real sigma2_def, real sigma2_min, |
808 | real *velectot, real *vvdwtot, real *dvdl) |
809 | { |
810 | real r, rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal; |
811 | real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2], sigma_powm2[2]; |
812 | real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw; |
813 | real rpinv, r_coul, r_vdw, velecsum, vvdwsum; |
814 | real fscal_vdw[2], fscal_elec[2]; |
815 | real velec[2], vvdw[2]; |
816 | int i, ntab; |
817 | |
818 | qq[0] = qqA; |
819 | qq[1] = qqB; |
820 | c6[0] = c6A; |
821 | c6[1] = c6B; |
822 | c12[0] = c12A; |
823 | c12[1] = c12B; |
824 | |
825 | if (sc_r_power == 6.0) |
826 | { |
827 | rpm2 = r2*r2; /* r4 */ |
828 | rp = rpm2*r2; /* r6 */ |
829 | } |
830 | else if (sc_r_power == 48.0) |
831 | { |
832 | rp = r2*r2*r2; /* r6 */ |
833 | rp = rp*rp; /* r12 */ |
834 | rp = rp*rp; /* r24 */ |
835 | rp = rp*rp; /* r48 */ |
836 | rpm2 = rp/r2; /* r46 */ |
837 | } |
838 | else |
839 | { |
840 | rp = pow(r2, 0.5*sc_r_power); /* not currently supported as input, but can handle it */ |
841 | rpm2 = rp/r2; |
842 | } |
843 | |
844 | /* Loop over state A(0) and B(1) */ |
845 | for (i = 0; i < 2; i++) |
846 | { |
847 | if ((c6[i] > 0) && (c12[i] > 0)) |
848 | { |
849 | /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively. |
850 | * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5. |
851 | */ |
852 | sigma6[i] = 0.5*c12[i]/c6[i]; |
853 | sigma2[i] = pow(0.5*c12[i]/c6[i], 1.0/3.0); |
854 | /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on |
855 | what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */ |
856 | if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */ |
857 | { |
858 | sigma6[i] = sigma6_min; |
859 | sigma2[i] = sigma2_min; |
860 | } |
861 | } |
862 | else |
863 | { |
864 | sigma6[i] = sigma6_def; |
865 | sigma2[i] = sigma2_def; |
866 | } |
867 | if (sc_r_power == 6.0) |
868 | { |
869 | sigma_pow[i] = sigma6[i]; |
870 | sigma_powm2[i] = sigma6[i]/sigma2[i]; |
871 | } |
872 | else if (sc_r_power == 48.0) |
873 | { |
874 | sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */ |
875 | sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */ |
876 | sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */ |
877 | sigma_powm2[i] = sigma_pow[i]/sigma2[i]; |
878 | } |
879 | else |
880 | { /* not really supported as input, but in here for testing the general case*/ |
881 | sigma_pow[i] = pow(sigma2[i], sc_r_power/2); |
882 | sigma_powm2[i] = sigma_pow[i]/(sigma2[i]); |
883 | } |
884 | } |
885 | |
886 | /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ |
887 | if ((c12[0] > 0) && (c12[1] > 0)) |
888 | { |
889 | alpha_vdw_eff = 0; |
890 | alpha_coul_eff = 0; |
891 | } |
892 | else |
893 | { |
894 | alpha_vdw_eff = alpha_vdw; |
895 | alpha_coul_eff = alpha_coul; |
896 | } |
897 | |
898 | /* Loop over A and B states again */ |
899 | for (i = 0; i < 2; i++) |
900 | { |
901 | fscal_elec[i] = 0; |
902 | fscal_vdw[i] = 0; |
903 | velec[i] = 0; |
904 | vvdw[i] = 0; |
905 | |
906 | /* Only spend time on A or B state if it is non-zero */ |
907 | if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) ) |
908 | { |
909 | /* Coulomb */ |
910 | rpinv = 1.0/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp); |
911 | r_coul = pow(rpinv, -1.0/sc_r_power); |
912 | |
913 | /* Electrostatics table lookup data */ |
914 | rtab = r_coul*tabscale; |
915 | ntab = rtab; |
916 | eps = rtab-ntab; |
917 | eps2 = eps*eps; |
918 | ntab = 12*ntab; |
919 | /* Electrostatics */ |
920 | Y = vftab[ntab]; |
921 | F = vftab[ntab+1]; |
922 | Geps = eps*vftab[ntab+2]; |
923 | Heps2 = eps2*vftab[ntab+3]; |
924 | Fp = F+Geps+Heps2; |
925 | VV = Y+eps*Fp; |
926 | FF = Fp+Geps+2.0*Heps2; |
927 | velec[i] = qq[i]*VV; |
928 | fscal_elec[i] = -qq[i]*FF*r_coul*rpinv*tabscale; |
929 | |
930 | /* Vdw */ |
931 | rpinv = 1.0/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp); |
932 | r_vdw = pow(rpinv, -1.0/sc_r_power); |
933 | /* Vdw table lookup data */ |
934 | rtab = r_vdw*tabscale; |
935 | ntab = rtab; |
936 | eps = rtab-ntab; |
937 | eps2 = eps*eps; |
938 | ntab = 12*ntab; |
939 | /* Dispersion */ |
940 | Y = vftab[ntab+4]; |
941 | F = vftab[ntab+5]; |
942 | Geps = eps*vftab[ntab+6]; |
943 | Heps2 = eps2*vftab[ntab+7]; |
944 | Fp = F+Geps+Heps2; |
945 | VV = Y+eps*Fp; |
946 | FF = Fp+Geps+2.0*Heps2; |
947 | vvdw[i] = c6[i]*VV; |
948 | fscal_vdw[i] = -c6[i]*FF; |
949 | |
950 | /* Repulsion */ |
951 | Y = vftab[ntab+8]; |
952 | F = vftab[ntab+9]; |
953 | Geps = eps*vftab[ntab+10]; |
954 | Heps2 = eps2*vftab[ntab+11]; |
955 | Fp = F+Geps+Heps2; |
956 | VV = Y+eps*Fp; |
957 | FF = Fp+Geps+2.0*Heps2; |
958 | vvdw[i] += c12[i]*VV; |
959 | fscal_vdw[i] -= c12[i]*FF; |
960 | fscal_vdw[i] *= r_vdw*rpinv*tabscale; |
961 | } |
962 | } |
963 | /* Now we have velec[i], vvdw[i], and fscal[i] for both states */ |
964 | /* Assemble A and B states */ |
965 | velecsum = 0; |
966 | vvdwsum = 0; |
967 | dvdl_coul = 0; |
968 | dvdl_vdw = 0; |
969 | fscal = 0; |
970 | for (i = 0; i < 2; i++) |
971 | { |
972 | velecsum += LFC[i]*velec[i]; |
973 | vvdwsum += LFV[i]*vvdw[i]; |
974 | |
975 | fscal += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2; |
976 | |
977 | dvdl_coul += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i]; |
978 | dvdl_vdw += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i]; |
979 | } |
980 | |
981 | dvdl[efptCOUL] += dvdl_coul; |
982 | dvdl[efptVDW] += dvdl_vdw; |
983 | |
984 | *velectot = velecsum; |
985 | *vvdwtot = vvdwsum; |
986 | |
987 | return fscal; |
988 | } |