File: | gromacs/mdlib/tables.c |
Location: | line 1435, column 5 |
Description: | Value stored to 'bGenTab' 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 | #include "gromacs/math/utilities.h" |
43 | #include "typedefs.h" |
44 | #include "names.h" |
45 | #include "gromacs/utility/smalloc.h" |
46 | #include "gromacs/utility/fatalerror.h" |
47 | #include "gromacs/utility/futil.h" |
48 | #include "gromacs/fileio/xvgr.h" |
49 | #include "gromacs/math/vec.h" |
50 | #include "network.h" |
51 | #include "physics.h" |
52 | #include "force.h" |
53 | #include "gromacs/fileio/gmxfio.h" |
54 | #include "macros.h" |
55 | #include "tables.h" |
56 | |
57 | /* All the possible (implemented) table functions */ |
58 | enum { |
59 | etabLJ6, |
60 | etabLJ12, |
61 | etabLJ6Shift, |
62 | etabLJ12Shift, |
63 | etabShift, |
64 | etabRF, |
65 | etabRF_ZERO, |
66 | etabCOUL, |
67 | etabEwald, |
68 | etabEwaldSwitch, |
69 | etabEwaldUser, |
70 | etabEwaldUserSwitch, |
71 | etabLJ6Ewald, |
72 | etabLJ6Switch, |
73 | etabLJ12Switch, |
74 | etabCOULSwitch, |
75 | etabLJ6Encad, |
76 | etabLJ12Encad, |
77 | etabCOULEncad, |
78 | etabEXPMIN, |
79 | etabUSER, |
80 | etabNR |
81 | }; |
82 | |
83 | /** Evaluates to true if the table type contains user data. */ |
84 | #define ETAB_USER(e)((e) == etabUSER || (e) == etabEwaldUser || (e) == etabEwaldUserSwitch ) ((e) == etabUSER || \ |
85 | (e) == etabEwaldUser || (e) == etabEwaldUserSwitch) |
86 | |
87 | typedef struct { |
88 | const char *name; |
89 | gmx_bool bCoulomb; |
90 | } t_tab_props; |
91 | |
92 | /* This structure holds name and a flag that tells whether |
93 | this is a Coulomb type funtion */ |
94 | static const t_tab_props tprops[etabNR] = { |
95 | { "LJ6", FALSE0 }, |
96 | { "LJ12", FALSE0 }, |
97 | { "LJ6Shift", FALSE0 }, |
98 | { "LJ12Shift", FALSE0 }, |
99 | { "Shift", TRUE1 }, |
100 | { "RF", TRUE1 }, |
101 | { "RF-zero", TRUE1 }, |
102 | { "COUL", TRUE1 }, |
103 | { "Ewald", TRUE1 }, |
104 | { "Ewald-Switch", TRUE1 }, |
105 | { "Ewald-User", TRUE1 }, |
106 | { "Ewald-User-Switch", TRUE1 }, |
107 | { "LJ6Ewald", FALSE0 }, |
108 | { "LJ6Switch", FALSE0 }, |
109 | { "LJ12Switch", FALSE0 }, |
110 | { "COULSwitch", TRUE1 }, |
111 | { "LJ6-Encad shift", FALSE0 }, |
112 | { "LJ12-Encad shift", FALSE0 }, |
113 | { "COUL-Encad shift", TRUE1 }, |
114 | { "EXPMIN", FALSE0 }, |
115 | { "USER", FALSE0 }, |
116 | }; |
117 | |
118 | /* Index in the table that says which function to use */ |
119 | enum { |
120 | etiCOUL, etiLJ6, etiLJ12, etiNR |
121 | }; |
122 | |
123 | typedef struct { |
124 | int nx, nx0; |
125 | double tabscale; |
126 | double *x, *v, *f; |
127 | } t_tabledata; |
128 | |
129 | #define pow2(x)((x)*(x)) ((x)*(x)) |
130 | #define pow3(x)((x)*(x)*(x)) ((x)*(x)*(x)) |
131 | #define pow4(x)((x)*(x)*(x)*(x)) ((x)*(x)*(x)*(x)) |
132 | #define pow5(x)((x)*(x)*(x)*(x)*(x)) ((x)*(x)*(x)*(x)*(x)) |
133 | |
134 | double v_q_ewald_lr(double beta, double r) |
135 | { |
136 | if (r == 0) |
137 | { |
138 | return beta*2/sqrt(M_PI3.14159265358979323846); |
139 | } |
140 | else |
141 | { |
142 | return gmx_erfd(beta*r)/r; |
143 | } |
144 | } |
145 | |
146 | double v_lj_ewald_lr(double beta, double r) |
147 | { |
148 | double br, br2, br4, r6, factor; |
149 | if (r == 0) |
150 | { |
151 | return pow(beta, 6)/6; |
152 | } |
153 | else |
154 | { |
155 | br = beta*r; |
156 | br2 = br*br; |
157 | br4 = br2*br2; |
158 | r6 = pow(r, 6.0); |
159 | factor = (1.0 - exp(-br2)*(1 + br2 + 0.5*br4))/r6; |
160 | return factor; |
161 | } |
162 | } |
163 | |
164 | void table_spline3_fill_ewald_lr(real *table_f, |
165 | real *table_v, |
166 | real *table_fdv0, |
167 | int ntab, |
168 | real dx, |
169 | real beta, |
170 | real_space_grid_contribution_computer v_lr) |
171 | { |
172 | real tab_max; |
173 | int i, i_inrange; |
174 | double dc, dc_new; |
175 | gmx_bool bOutOfRange; |
176 | double v_r0, v_r1, v_inrange, vi, a0, a1, a2dx; |
177 | double x_r0; |
178 | |
179 | /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument |
180 | * depending on wether we should create electrostatic or Lennard-Jones Ewald tables. |
181 | */ |
182 | |
183 | if (ntab < 2) |
184 | { |
185 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 185, "Can not make a spline table with less than 2 points"); |
186 | } |
187 | |
188 | /* We need some margin to be able to divide table values by r |
189 | * in the kernel and also to do the integration arithmetics |
190 | * without going out of range. Furthemore, we divide by dx below. |
191 | */ |
192 | tab_max = GMX_REAL_MAX3.40282346E+38*0.0001; |
193 | |
194 | /* This function produces a table with: |
195 | * maximum energy error: V'''/(6*12*sqrt(3))*dx^3 |
196 | * maximum force error: V'''/(6*4)*dx^2 |
197 | * The rms force error is the max error times 1/sqrt(5)=0.45. |
198 | */ |
199 | |
200 | bOutOfRange = FALSE0; |
201 | i_inrange = ntab; |
202 | v_inrange = 0; |
203 | dc = 0; |
204 | for (i = ntab-1; i >= 0; i--) |
205 | { |
206 | x_r0 = i*dx; |
207 | |
208 | v_r0 = (*v_lr)(beta, x_r0); |
209 | |
210 | if (!bOutOfRange) |
211 | { |
212 | i_inrange = i; |
213 | v_inrange = v_r0; |
214 | |
215 | vi = v_r0; |
216 | } |
217 | else |
218 | { |
219 | /* Linear continuation for the last point in range */ |
220 | vi = v_inrange - dc*(i - i_inrange)*dx; |
221 | } |
222 | |
223 | if (table_v != NULL((void*)0)) |
224 | { |
225 | table_v[i] = vi; |
226 | } |
227 | |
228 | if (i == 0) |
229 | { |
230 | continue; |
231 | } |
232 | |
233 | /* Get the potential at table point i-1 */ |
234 | v_r1 = (*v_lr)(beta, (i-1)*dx); |
235 | |
236 | if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max) |
237 | { |
238 | bOutOfRange = TRUE1; |
239 | } |
240 | |
241 | if (!bOutOfRange) |
242 | { |
243 | /* Calculate the average second derivative times dx over interval i-1 to i. |
244 | * Using the function values at the end points and in the middle. |
245 | */ |
246 | a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx); |
247 | /* Set the derivative of the spline to match the difference in potential |
248 | * over the interval plus the average effect of the quadratic term. |
249 | * This is the essential step for minimizing the error in the force. |
250 | */ |
251 | dc = (v_r0 - v_r1)/dx + 0.5*a2dx; |
252 | } |
253 | |
254 | if (i == ntab - 1) |
255 | { |
256 | /* Fill the table with the force, minus the derivative of the spline */ |
257 | table_f[i] = -dc; |
258 | } |
259 | else |
260 | { |
261 | /* tab[i] will contain the average of the splines over the two intervals */ |
262 | table_f[i] += -0.5*dc; |
263 | } |
264 | |
265 | if (!bOutOfRange) |
266 | { |
267 | /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2 |
268 | * matching the potential at the two end points |
269 | * and the derivative dc at the end point xr. |
270 | */ |
271 | a0 = v_r0; |
272 | a1 = dc; |
273 | a2dx = (a1*dx + v_r1 - a0)*2/dx; |
274 | |
275 | /* Set dc to the derivative at the next point */ |
276 | dc_new = a1 - a2dx; |
277 | |
278 | if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max) |
279 | { |
280 | bOutOfRange = TRUE1; |
281 | } |
282 | else |
283 | { |
284 | dc = dc_new; |
285 | } |
286 | } |
287 | |
288 | table_f[(i-1)] = -0.5*dc; |
289 | } |
290 | /* Currently the last value only contains half the force: double it */ |
291 | table_f[0] *= 2; |
292 | |
293 | if (table_v != NULL((void*)0) && table_fdv0 != NULL((void*)0)) |
294 | { |
295 | /* Copy to FDV0 table too. Allocation occurs in forcerec.c, |
296 | * init_ewald_f_table(). |
297 | */ |
298 | for (i = 0; i < ntab-1; i++) |
299 | { |
300 | table_fdv0[4*i] = table_f[i]; |
301 | table_fdv0[4*i+1] = table_f[i+1]-table_f[i]; |
302 | table_fdv0[4*i+2] = table_v[i]; |
303 | table_fdv0[4*i+3] = 0.0; |
304 | } |
305 | table_fdv0[4*(ntab-1)] = table_f[(ntab-1)]; |
306 | table_fdv0[4*(ntab-1)+1] = -table_f[(ntab-1)]; |
307 | table_fdv0[4*(ntab-1)+2] = table_v[(ntab-1)]; |
308 | table_fdv0[4*(ntab-1)+3] = 0.0; |
309 | } |
310 | } |
311 | |
312 | /* The scale (1/spacing) for third order spline interpolation |
313 | * of the Ewald mesh contribution which needs to be subtracted |
314 | * from the non-bonded interactions. |
315 | */ |
316 | real ewald_spline3_table_scale(real ewaldcoeff, real rc) |
317 | { |
318 | double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */ |
319 | double ftol, etol; |
320 | double sc_f, sc_e; |
321 | |
322 | /* Force tolerance: single precision accuracy */ |
323 | ftol = GMX_FLOAT_EPS5.96046448E-08; |
324 | sc_f = sqrt(erf_x_d3/(6*4*ftol*ewaldcoeff))*ewaldcoeff; |
325 | |
326 | /* Energy tolerance: 10x more accurate than the cut-off jump */ |
327 | etol = 0.1*gmx_erfc(ewaldcoeff*rc)gmx_erfcf(ewaldcoeff*rc); |
328 | etol = max(etol, GMX_REAL_EPS)(((etol) > (5.96046448E-08)) ? (etol) : (5.96046448E-08) ); |
329 | sc_e = pow(erf_x_d3/(6*12*sqrt(3)*etol), 1.0/3.0)*ewaldcoeff; |
330 | |
331 | return max(sc_f, sc_e)(((sc_f) > (sc_e)) ? (sc_f) : (sc_e) ); |
332 | } |
333 | |
334 | /* Calculate the potential and force for an r value |
335 | * in exactly the same way it is done in the inner loop. |
336 | * VFtab is a pointer to the table data, offset is |
337 | * the point where we should begin and stride is |
338 | * 4 if we have a buckingham table, 3 otherwise. |
339 | * If you want to evaluate table no N, set offset to 4*N. |
340 | * |
341 | * We use normal precision here, since that is what we |
342 | * will use in the inner loops. |
343 | */ |
344 | static void evaluate_table(real VFtab[], int offset, int stride, |
345 | real tabscale, real r, real *y, real *yp) |
346 | { |
347 | int n; |
348 | real rt, eps, eps2; |
349 | real Y, F, Geps, Heps2, Fp; |
350 | |
351 | rt = r*tabscale; |
352 | n = (int)rt; |
353 | eps = rt - n; |
354 | eps2 = eps*eps; |
355 | n = offset+stride*n; |
356 | Y = VFtab[n]; |
357 | F = VFtab[n+1]; |
358 | Geps = eps*VFtab[n+2]; |
359 | Heps2 = eps2*VFtab[n+3]; |
360 | Fp = F+Geps+Heps2; |
361 | *y = Y+eps*Fp; |
362 | *yp = (Fp+Geps+2.0*Heps2)*tabscale; |
363 | } |
364 | |
365 | static void copy2table(int n, int offset, int stride, |
366 | double x[], double Vtab[], double Ftab[], real scalefactor, |
367 | real dest[]) |
368 | { |
369 | /* Use double prec. for the intermediary variables |
370 | * and temporary x/vtab/vtab2 data to avoid unnecessary |
371 | * loss of precision. |
372 | */ |
373 | int i, nn0; |
374 | double F, G, H, h; |
375 | |
376 | h = 0; |
377 | for (i = 0; (i < n); i++) |
378 | { |
379 | if (i < n-1) |
380 | { |
381 | h = x[i+1] - x[i]; |
382 | F = -Ftab[i]*h; |
383 | G = 3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h; |
384 | H = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] + Ftab[i])*h; |
385 | } |
386 | else |
387 | { |
388 | /* Fill the last entry with a linear potential, |
389 | * this is mainly for rounding issues with angle and dihedral potentials. |
390 | */ |
391 | F = -Ftab[i]*h; |
392 | G = 0; |
393 | H = 0; |
394 | } |
395 | nn0 = offset + i*stride; |
396 | dest[nn0] = scalefactor*Vtab[i]; |
397 | dest[nn0+1] = scalefactor*F; |
398 | dest[nn0+2] = scalefactor*G; |
399 | dest[nn0+3] = scalefactor*H; |
400 | } |
401 | } |
402 | |
403 | static void init_table(int n, int nx0, |
404 | double tabscale, t_tabledata *td, gmx_bool bAlloc) |
405 | { |
406 | int i; |
407 | |
408 | td->nx = n; |
409 | td->nx0 = nx0; |
410 | td->tabscale = tabscale; |
411 | if (bAlloc) |
412 | { |
413 | snew(td->x, td->nx)(td->x) = save_calloc("td->x", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 413, (td->nx), sizeof(*(td->x))); |
414 | snew(td->v, td->nx)(td->v) = save_calloc("td->v", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 414, (td->nx), sizeof(*(td->v))); |
415 | snew(td->f, td->nx)(td->f) = save_calloc("td->f", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 415, (td->nx), sizeof(*(td->f))); |
416 | } |
417 | for (i = 0; (i < td->nx); i++) |
418 | { |
419 | td->x[i] = i/tabscale; |
420 | } |
421 | } |
422 | |
423 | static void spline_forces(int nx, double h, double v[], gmx_bool bS3, gmx_bool bE3, |
424 | double f[]) |
425 | { |
426 | int start, end, i; |
427 | double v3, b_s, b_e, b; |
428 | double beta, *gamma; |
429 | |
430 | /* Formulas can be found in: |
431 | * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007 |
432 | */ |
433 | |
434 | if (nx < 4 && (bS3 || bE3)) |
435 | { |
436 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 436, "Can not generate splines with third derivative boundary conditions with less than 4 (%d) points", nx); |
437 | } |
438 | |
439 | /* To make life easy we initially set the spacing to 1 |
440 | * and correct for this at the end. |
441 | */ |
442 | beta = 2; |
443 | if (bS3) |
444 | { |
445 | /* Fit V''' at the start */ |
446 | v3 = v[3] - 3*v[2] + 3*v[1] - v[0]; |
447 | if (debug) |
448 | { |
449 | fprintf(debug, "The left third derivative is %g\n", v3/(h*h*h)); |
450 | } |
451 | b_s = 2*(v[1] - v[0]) + v3/6; |
452 | start = 0; |
453 | |
454 | if (FALSE0) |
455 | { |
456 | /* Fit V'' at the start */ |
457 | real v2; |
458 | |
459 | v2 = -v[3] + 4*v[2] - 5*v[1] + 2*v[0]; |
460 | /* v2 = v[2] - 2*v[1] + v[0]; */ |
461 | if (debug) |
462 | { |
463 | fprintf(debug, "The left second derivative is %g\n", v2/(h*h)); |
464 | } |
465 | b_s = 3*(v[1] - v[0]) - v2/2; |
466 | start = 0; |
467 | } |
468 | } |
469 | else |
470 | { |
471 | b_s = 3*(v[2] - v[0]) + f[0]*h; |
472 | start = 1; |
473 | } |
474 | if (bE3) |
475 | { |
476 | /* Fit V''' at the end */ |
477 | v3 = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4]; |
478 | if (debug) |
479 | { |
480 | fprintf(debug, "The right third derivative is %g\n", v3/(h*h*h)); |
481 | } |
482 | b_e = 2*(v[nx-1] - v[nx-2]) + v3/6; |
483 | end = nx; |
484 | } |
485 | else |
486 | { |
487 | /* V'=0 at the end */ |
488 | b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h; |
489 | end = nx - 1; |
490 | } |
491 | |
492 | snew(gamma, nx)(gamma) = save_calloc("gamma", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 492, (nx), sizeof(*(gamma))); |
493 | beta = (bS3 ? 1 : 4); |
494 | |
495 | /* For V'' fitting */ |
496 | /* beta = (bS3 ? 2 : 4); */ |
497 | |
498 | f[start] = b_s/beta; |
499 | for (i = start+1; i < end; i++) |
500 | { |
501 | gamma[i] = 1/beta; |
502 | beta = 4 - gamma[i]; |
503 | b = 3*(v[i+1] - v[i-1]); |
504 | f[i] = (b - f[i-1])/beta; |
505 | } |
506 | gamma[end-1] = 1/beta; |
507 | beta = (bE3 ? 1 : 4) - gamma[end-1]; |
508 | f[end-1] = (b_e - f[end-2])/beta; |
509 | |
510 | for (i = end-2; i >= start; i--) |
511 | { |
512 | f[i] -= gamma[i+1]*f[i+1]; |
513 | } |
514 | sfree(gamma)save_free("gamma", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 514, (gamma)); |
515 | |
516 | /* Correct for the minus sign and the spacing */ |
517 | for (i = start; i < end; i++) |
518 | { |
519 | f[i] = -f[i]/h; |
520 | } |
521 | } |
522 | |
523 | static void set_forces(FILE *fp, int angle, |
524 | int nx, double h, double v[], double f[], |
525 | int table) |
526 | { |
527 | int start, end; |
528 | |
529 | if (angle == 2) |
530 | { |
531 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 531, |
532 | "Force generation for dihedral tables is not (yet) implemented"); |
533 | } |
534 | |
535 | start = 0; |
536 | while (v[start] == 0) |
537 | { |
538 | start++; |
539 | } |
540 | |
541 | end = nx; |
542 | while (v[end-1] == 0) |
543 | { |
544 | end--; |
545 | } |
546 | if (end > nx - 2) |
547 | { |
548 | end = nx; |
549 | } |
550 | else |
551 | { |
552 | end++; |
553 | } |
554 | |
555 | if (fp) |
556 | { |
557 | fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n", |
558 | table+1, start*h, end == nx ? "V'''" : "V'=0", (end-1)*h); |
559 | } |
560 | spline_forces(end-start, h, v+start, TRUE1, end == nx, f+start); |
561 | } |
562 | |
563 | static void read_tables(FILE *fp, const char *fn, |
564 | int ntab, int angle, t_tabledata td[]) |
565 | { |
566 | char *libfn; |
567 | char buf[STRLEN4096]; |
568 | double **yy = NULL((void*)0), start, end, dx0, dx1, ssd, vm, vp, f, numf; |
569 | int k, i, nx, nx0 = 0, ny, nny, ns; |
570 | gmx_bool bAllZero, bZeroV, bZeroF; |
571 | double tabscale; |
572 | |
573 | nny = 2*ntab+1; |
574 | libfn = gmxlibfn(fn); |
575 | nx = read_xvg(libfn, &yy, &ny); |
576 | if (ny != nny) |
577 | { |
578 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 578, "Trying to read file %s, but nr columns = %d, should be %d", |
579 | libfn, ny, nny); |
580 | } |
581 | if (angle == 0) |
582 | { |
583 | if (yy[0][0] != 0.0) |
584 | { |
585 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 585, |
586 | "The first distance in file %s is %f nm instead of %f nm", |
587 | libfn, yy[0][0], 0.0); |
588 | } |
589 | } |
590 | else |
591 | { |
592 | if (angle == 1) |
593 | { |
594 | start = 0.0; |
595 | } |
596 | else |
597 | { |
598 | start = -180.0; |
599 | } |
600 | end = 180.0; |
601 | if (yy[0][0] != start || yy[0][nx-1] != end) |
602 | { |
603 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 603, "The angles in file %s should go from %f to %f instead of %f to %f\n", |
604 | libfn, start, end, yy[0][0], yy[0][nx-1]); |
605 | } |
606 | } |
607 | |
608 | tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]); |
609 | |
610 | if (fp) |
611 | { |
612 | fprintf(fp, "Read user tables from %s with %d data points.\n", libfn, nx); |
613 | if (angle == 0) |
614 | { |
615 | fprintf(fp, "Tabscale = %g points/nm\n", tabscale); |
616 | } |
617 | } |
618 | |
619 | bAllZero = TRUE1; |
620 | for (k = 0; k < ntab; k++) |
621 | { |
622 | bZeroV = TRUE1; |
623 | bZeroF = TRUE1; |
624 | for (i = 0; (i < nx); i++) |
625 | { |
626 | if (i >= 2) |
627 | { |
628 | dx0 = yy[0][i-1] - yy[0][i-2]; |
629 | dx1 = yy[0][i] - yy[0][i-1]; |
630 | /* Check for 1% deviation in spacing */ |
631 | if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1))) |
632 | { |
633 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 633, "In table file '%s' the x values are not equally spaced: %f %f %f", fn, yy[0][i-2], yy[0][i-1], yy[0][i]); |
634 | } |
635 | } |
636 | if (yy[1+k*2][i] != 0) |
637 | { |
638 | bZeroV = FALSE0; |
639 | if (bAllZero) |
640 | { |
641 | bAllZero = FALSE0; |
642 | nx0 = i; |
643 | } |
644 | if (yy[1+k*2][i] > 0.01*GMX_REAL_MAX3.40282346E+38 || |
645 | yy[1+k*2][i] < -0.01*GMX_REAL_MAX3.40282346E+38) |
646 | { |
647 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 647, "Out of range potential value %g in file '%s'", |
648 | yy[1+k*2][i], fn); |
649 | } |
650 | } |
651 | if (yy[1+k*2+1][i] != 0) |
652 | { |
653 | bZeroF = FALSE0; |
654 | if (bAllZero) |
655 | { |
656 | bAllZero = FALSE0; |
657 | nx0 = i; |
658 | } |
659 | if (yy[1+k*2+1][i] > 0.01*GMX_REAL_MAX3.40282346E+38 || |
660 | yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX3.40282346E+38) |
661 | { |
662 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 662, "Out of range force value %g in file '%s'", |
663 | yy[1+k*2+1][i], fn); |
664 | } |
665 | } |
666 | } |
667 | |
668 | if (!bZeroV && bZeroF) |
669 | { |
670 | set_forces(fp, angle, nx, 1/tabscale, yy[1+k*2], yy[1+k*2+1], k); |
671 | } |
672 | else |
673 | { |
674 | /* Check if the second column is close to minus the numerical |
675 | * derivative of the first column. |
676 | */ |
677 | ssd = 0; |
678 | ns = 0; |
679 | for (i = 1; (i < nx-1); i++) |
680 | { |
681 | vm = yy[1+2*k][i-1]; |
682 | vp = yy[1+2*k][i+1]; |
683 | f = yy[1+2*k+1][i]; |
684 | if (vm != 0 && vp != 0 && f != 0) |
685 | { |
686 | /* Take the centered difference */ |
687 | numf = -(vp - vm)*0.5*tabscale; |
688 | ssd += fabs(2*(f - numf)/(f + numf)); |
689 | ns++; |
690 | } |
691 | } |
692 | if (ns > 0) |
693 | { |
694 | ssd /= ns; |
695 | sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %d%% from minus the numerical derivative of the potential\n", ns, k, libfn, (int)(100*ssd+0.5)); |
696 | if (debug) |
697 | { |
698 | fprintf(debug, "%s", buf); |
699 | } |
700 | if (ssd > 0.2) |
701 | { |
702 | if (fp) |
703 | { |
704 | fprintf(fp, "\nWARNING: %s\n", buf); |
705 | } |
706 | fprintf(stderrstderr, "\nWARNING: %s\n", buf); |
707 | } |
708 | } |
709 | } |
710 | } |
711 | if (bAllZero && fp) |
712 | { |
713 | fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn); |
714 | } |
715 | |
716 | for (k = 0; (k < ntab); k++) |
717 | { |
718 | init_table(nx, nx0, tabscale, &(td[k]), TRUE1); |
719 | for (i = 0; (i < nx); i++) |
720 | { |
721 | td[k].x[i] = yy[0][i]; |
722 | td[k].v[i] = yy[2*k+1][i]; |
723 | td[k].f[i] = yy[2*k+2][i]; |
724 | } |
725 | } |
726 | for (i = 0; (i < ny); i++) |
727 | { |
728 | sfree(yy[i])save_free("yy[i]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 728, (yy[i])); |
729 | } |
730 | sfree(yy)save_free("yy", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 730, (yy)); |
731 | sfree(libfn)save_free("libfn", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 731, (libfn)); |
732 | } |
733 | |
734 | static void done_tabledata(t_tabledata *td) |
735 | { |
736 | int i; |
737 | |
738 | if (!td) |
739 | { |
740 | return; |
741 | } |
742 | |
743 | sfree(td->x)save_free("td->x", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 743, (td->x)); |
744 | sfree(td->v)save_free("td->v", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 744, (td->v)); |
745 | sfree(td->f)save_free("td->f", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 745, (td->f)); |
746 | } |
747 | |
748 | static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr) |
749 | { |
750 | /* Fill the table according to the formulas in the manual. |
751 | * In principle, we only need the potential and the second |
752 | * derivative, but then we would have to do lots of calculations |
753 | * in the inner loop. By precalculating some terms (see manual) |
754 | * we get better eventual performance, despite a larger table. |
755 | * |
756 | * Since some of these higher-order terms are very small, |
757 | * we always use double precision to calculate them here, in order |
758 | * to avoid unnecessary loss of precision. |
759 | */ |
760 | #ifdef DEBUG_SWITCH |
761 | FILE *fp; |
762 | #endif |
763 | int i; |
764 | double reppow, p; |
765 | double r1, rc, r12, r13; |
766 | double r, r2, r6, rc6; |
767 | double expr, Vtab, Ftab; |
768 | /* Parameters for David's function */ |
769 | double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0; |
770 | /* Parameters for the switching function */ |
771 | double ksw, swi, swi1; |
772 | /* Temporary parameters */ |
773 | gmx_bool bSwitch, bShift; |
774 | double ewc = fr->ewaldcoeff_q; |
775 | double ewclj = fr->ewaldcoeff_lj; |
776 | |
777 | bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) || |
778 | (tp == etabCOULSwitch) || |
779 | (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch)); |
780 | |
781 | bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) || |
782 | (tp == etabShift)); |
783 | |
784 | reppow = fr->reppow; |
785 | |
786 | if (tprops[tp].bCoulomb) |
787 | { |
788 | r1 = fr->rcoulomb_switch; |
789 | rc = fr->rcoulomb; |
790 | } |
791 | else |
792 | { |
793 | r1 = fr->rvdw_switch; |
794 | rc = fr->rvdw; |
795 | } |
796 | if (bSwitch) |
797 | { |
798 | ksw = 1.0/(pow5(rc-r1)((rc-r1)*(rc-r1)*(rc-r1)*(rc-r1)*(rc-r1))); |
799 | } |
800 | else |
801 | { |
802 | ksw = 0.0; |
803 | } |
804 | if (bShift) |
805 | { |
806 | if (tp == etabShift) |
807 | { |
808 | p = 1; |
809 | } |
810 | else if (tp == etabLJ6Shift) |
811 | { |
812 | p = 6; |
813 | } |
814 | else |
815 | { |
816 | p = reppow; |
817 | } |
818 | |
819 | A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc, p+2)*pow2(rc-r1)((rc-r1)*(rc-r1))); |
820 | B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc, p+2)*pow3(rc-r1)((rc-r1)*(rc-r1)*(rc-r1))); |
821 | C = 1.0/pow(rc, p)-A/3.0*pow3(rc-r1)((rc-r1)*(rc-r1)*(rc-r1))-B/4.0*pow4(rc-r1)((rc-r1)*(rc-r1)*(rc-r1)*(rc-r1)); |
822 | if (tp == etabLJ6Shift) |
823 | { |
824 | A = -A; |
825 | B = -B; |
826 | C = -C; |
827 | } |
828 | A_3 = A/3.0; |
829 | B_4 = B/4.0; |
830 | } |
831 | if (debug) |
832 | { |
833 | fprintf(debug, "Setting up tables\n"); fflush(debug); |
834 | } |
835 | |
836 | #ifdef DEBUG_SWITCH |
837 | fp = xvgropen("switch.xvg", "switch", "r", "s"); |
838 | #endif |
839 | |
840 | for (i = td->nx0; (i < td->nx); i++) |
841 | { |
842 | r = td->x[i]; |
843 | r2 = r*r; |
844 | r6 = 1.0/(r2*r2*r2); |
845 | if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS1.11022302E-16)) |
846 | { |
847 | r12 = r6*r6; |
848 | } |
849 | else |
850 | { |
851 | r12 = pow(r, -reppow); |
852 | } |
853 | Vtab = 0.0; |
854 | Ftab = 0.0; |
855 | if (bSwitch) |
856 | { |
857 | /* swi is function, swi1 1st derivative and swi2 2nd derivative */ |
858 | /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for |
859 | * r1<=r<=rc. The 1st and 2nd derivatives are both zero at |
860 | * r1 and rc. |
861 | * ksw is just the constant 1/(rc-r1)^5, to save some calculations... |
862 | */ |
863 | if (r <= r1) |
864 | { |
865 | swi = 1.0; |
866 | swi1 = 0.0; |
867 | } |
868 | else if (r >= rc) |
869 | { |
870 | swi = 0.0; |
871 | swi1 = 0.0; |
872 | } |
873 | else |
874 | { |
875 | swi = 1 - 10*pow3(r-r1)((r-r1)*(r-r1)*(r-r1))*ksw*pow2(rc-r1)((rc-r1)*(rc-r1)) |
876 | + 15*pow4(r-r1)((r-r1)*(r-r1)*(r-r1)*(r-r1))*ksw*(rc-r1) - 6*pow5(r-r1)((r-r1)*(r-r1)*(r-r1)*(r-r1)*(r-r1))*ksw; |
877 | swi1 = -30*pow2(r-r1)((r-r1)*(r-r1))*ksw*pow2(rc-r1)((rc-r1)*(rc-r1)) |
878 | + 60*pow3(r-r1)((r-r1)*(r-r1)*(r-r1))*ksw*(rc-r1) - 30*pow4(r-r1)((r-r1)*(r-r1)*(r-r1)*(r-r1))*ksw; |
879 | } |
880 | } |
881 | else /* not really needed, but avoids compiler warnings... */ |
882 | { |
883 | swi = 1.0; |
884 | swi1 = 0.0; |
885 | } |
886 | #ifdef DEBUG_SWITCH |
887 | fprintf(fp, "%10g %10g %10g %10g\n", r, swi, swi1, swi2); |
888 | #endif |
889 | |
890 | rc6 = rc*rc*rc; |
891 | rc6 = 1.0/(rc6*rc6); |
892 | |
893 | switch (tp) |
894 | { |
895 | case etabLJ6: |
896 | /* Dispersion */ |
897 | Vtab = -r6; |
898 | Ftab = 6.0*Vtab/r; |
899 | break; |
900 | case etabLJ6Switch: |
901 | case etabLJ6Shift: |
902 | /* Dispersion */ |
903 | if (r < rc) |
904 | { |
905 | Vtab = -r6; |
906 | Ftab = 6.0*Vtab/r; |
907 | break; |
908 | } |
909 | break; |
910 | case etabLJ12: |
911 | /* Repulsion */ |
912 | Vtab = r12; |
913 | Ftab = reppow*Vtab/r; |
914 | break; |
915 | case etabLJ12Switch: |
916 | case etabLJ12Shift: |
917 | /* Repulsion */ |
918 | if (r < rc) |
919 | { |
920 | Vtab = r12; |
921 | Ftab = reppow*Vtab/r; |
922 | } |
923 | break; |
924 | case etabLJ6Encad: |
925 | if (r < rc) |
926 | { |
927 | Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6); |
928 | Ftab = -(6.0*r6/r-6.0*rc6/rc); |
929 | } |
930 | else /* r>rc */ |
931 | { |
932 | Vtab = 0; |
933 | Ftab = 0; |
934 | } |
935 | break; |
936 | case etabLJ12Encad: |
937 | if (r < rc) |
938 | { |
939 | Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6); |
940 | Ftab = -(6.0*r6/r-6.0*rc6/rc); |
941 | } |
942 | else /* r>rc */ |
943 | { |
944 | Vtab = 0; |
945 | Ftab = 0; |
946 | } |
947 | break; |
948 | case etabCOUL: |
949 | Vtab = 1.0/r; |
950 | Ftab = 1.0/r2; |
951 | break; |
952 | case etabCOULSwitch: |
953 | case etabShift: |
954 | if (r < rc) |
955 | { |
956 | Vtab = 1.0/r; |
957 | Ftab = 1.0/r2; |
958 | } |
959 | break; |
960 | case etabEwald: |
961 | case etabEwaldSwitch: |
962 | Vtab = gmx_erfc(ewc*r)gmx_erfcf(ewc*r)/r; |
963 | Ftab = gmx_erfc(ewc*r)gmx_erfcf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI1.12837916709551257390/r; |
964 | break; |
965 | case etabEwaldUser: |
966 | case etabEwaldUserSwitch: |
967 | /* Only calculate the negative of the reciprocal space contribution */ |
968 | Vtab = -gmx_erf(ewc*r)gmx_erff(ewc*r)/r; |
969 | Ftab = -gmx_erf(ewc*r)gmx_erff(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI1.12837916709551257390/r; |
970 | break; |
971 | case etabLJ6Ewald: |
972 | Vtab = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + pow4(ewclj)((ewclj)*(ewclj)*(ewclj)*(ewclj))*r2*r2/2); |
973 | Ftab = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*pow5(ewclj)((ewclj)*(ewclj)*(ewclj)*(ewclj)*(ewclj))*ewclj*r2*r2*r; |
974 | break; |
975 | case etabRF: |
976 | case etabRF_ZERO: |
977 | Vtab = 1.0/r + fr->k_rf*r2 - fr->c_rf; |
978 | Ftab = 1.0/r2 - 2*fr->k_rf*r; |
979 | if (tp == etabRF_ZERO && r >= rc) |
980 | { |
981 | Vtab = 0; |
982 | Ftab = 0; |
983 | } |
984 | break; |
985 | case etabEXPMIN: |
986 | expr = exp(-r); |
987 | Vtab = expr; |
988 | Ftab = expr; |
989 | break; |
990 | case etabCOULEncad: |
991 | if (r < rc) |
992 | { |
993 | Vtab = 1.0/r-(rc-r)/(rc*rc)-1.0/rc; |
994 | Ftab = 1.0/r2-1.0/(rc*rc); |
995 | } |
996 | else /* r>rc */ |
997 | { |
998 | Vtab = 0; |
999 | Ftab = 0; |
1000 | } |
1001 | break; |
1002 | default: |
1003 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1003, "Table type %d not implemented yet. (%s,%d)", |
1004 | tp, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", __LINE__1004); |
1005 | } |
1006 | if (bShift) |
1007 | { |
1008 | /* Normal coulomb with cut-off correction for potential */ |
1009 | if (r < rc) |
1010 | { |
1011 | Vtab -= C; |
1012 | /* If in Shifting range add something to it */ |
1013 | if (r > r1) |
1014 | { |
1015 | r12 = (r-r1)*(r-r1); |
1016 | r13 = (r-r1)*r12; |
1017 | Vtab += -A_3*r13 - B_4*r12*r12; |
1018 | Ftab += A*r12 + B*r13; |
1019 | } |
1020 | } |
1021 | } |
1022 | |
1023 | if (ETAB_USER(tp)((tp) == etabUSER || (tp) == etabEwaldUser || (tp) == etabEwaldUserSwitch )) |
1024 | { |
1025 | Vtab += td->v[i]; |
1026 | Ftab += td->f[i]; |
1027 | } |
1028 | |
1029 | if ((r > r1) && bSwitch) |
1030 | { |
1031 | Ftab = Ftab*swi - Vtab*swi1; |
1032 | Vtab = Vtab*swi; |
1033 | } |
1034 | |
1035 | /* Convert to single precision when we store to mem */ |
1036 | td->v[i] = Vtab; |
1037 | td->f[i] = Ftab; |
1038 | } |
1039 | |
1040 | /* Continue the table linearly from nx0 to 0. |
1041 | * These values are only required for energy minimization with overlap or TPI. |
1042 | */ |
1043 | for (i = td->nx0-1; i >= 0; i--) |
1044 | { |
1045 | td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]); |
1046 | td->f[i] = td->f[i+1]; |
1047 | } |
1048 | |
1049 | #ifdef DEBUG_SWITCH |
1050 | gmx_fio_fclose(fp); |
1051 | #endif |
1052 | } |
1053 | |
1054 | static void set_table_type(int tabsel[], const t_forcerec *fr, gmx_bool b14only) |
1055 | { |
1056 | int eltype, vdwtype; |
1057 | |
1058 | /* Set the different table indices. |
1059 | * Coulomb first. |
1060 | */ |
1061 | |
1062 | |
1063 | if (b14only) |
1064 | { |
1065 | switch (fr->eeltype) |
1066 | { |
1067 | case eelRF_NEC: |
1068 | eltype = eelRF; |
1069 | break; |
1070 | case eelUSER: |
1071 | case eelPMEUSER: |
1072 | case eelPMEUSERSWITCH: |
1073 | eltype = eelUSER; |
1074 | break; |
1075 | default: |
1076 | eltype = eelCUT; |
1077 | } |
1078 | } |
1079 | else |
1080 | { |
1081 | eltype = fr->eeltype; |
1082 | } |
1083 | |
1084 | switch (eltype) |
1085 | { |
1086 | case eelCUT: |
1087 | tabsel[etiCOUL] = etabCOUL; |
1088 | break; |
1089 | case eelPOISSON: |
1090 | tabsel[etiCOUL] = etabShift; |
1091 | break; |
1092 | case eelSHIFT: |
1093 | if (fr->rcoulomb > fr->rcoulomb_switch) |
1094 | { |
1095 | tabsel[etiCOUL] = etabShift; |
1096 | } |
1097 | else |
1098 | { |
1099 | tabsel[etiCOUL] = etabCOUL; |
1100 | } |
1101 | break; |
1102 | case eelEWALD: |
1103 | case eelPME: |
1104 | case eelP3M_AD: |
1105 | tabsel[etiCOUL] = etabEwald; |
1106 | break; |
1107 | case eelPMESWITCH: |
1108 | tabsel[etiCOUL] = etabEwaldSwitch; |
1109 | break; |
1110 | case eelPMEUSER: |
1111 | tabsel[etiCOUL] = etabEwaldUser; |
1112 | break; |
1113 | case eelPMEUSERSWITCH: |
1114 | tabsel[etiCOUL] = etabEwaldUserSwitch; |
1115 | break; |
1116 | case eelRF: |
1117 | case eelGRF: |
1118 | case eelRF_NEC: |
1119 | tabsel[etiCOUL] = etabRF; |
1120 | break; |
1121 | case eelRF_ZERO: |
1122 | tabsel[etiCOUL] = etabRF_ZERO; |
1123 | break; |
1124 | case eelSWITCH: |
1125 | tabsel[etiCOUL] = etabCOULSwitch; |
1126 | break; |
1127 | case eelUSER: |
1128 | tabsel[etiCOUL] = etabUSER; |
1129 | break; |
1130 | case eelENCADSHIFT: |
1131 | tabsel[etiCOUL] = etabCOULEncad; |
1132 | break; |
1133 | default: |
1134 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1134, "Invalid eeltype %d", eltype); |
1135 | } |
1136 | |
1137 | /* Van der Waals time */ |
1138 | if (fr->bBHAM && !b14only) |
1139 | { |
1140 | tabsel[etiLJ6] = etabLJ6; |
1141 | tabsel[etiLJ12] = etabEXPMIN; |
1142 | } |
1143 | else |
1144 | { |
1145 | if (b14only && fr->vdwtype != evdwUSER) |
1146 | { |
1147 | vdwtype = evdwCUT; |
1148 | } |
1149 | else |
1150 | { |
1151 | vdwtype = fr->vdwtype; |
1152 | } |
1153 | |
1154 | switch (vdwtype) |
1155 | { |
1156 | case evdwSWITCH: |
1157 | tabsel[etiLJ6] = etabLJ6Switch; |
1158 | tabsel[etiLJ12] = etabLJ12Switch; |
1159 | break; |
1160 | case evdwSHIFT: |
1161 | tabsel[etiLJ6] = etabLJ6Shift; |
1162 | tabsel[etiLJ12] = etabLJ12Shift; |
1163 | break; |
1164 | case evdwUSER: |
1165 | tabsel[etiLJ6] = etabUSER; |
1166 | tabsel[etiLJ12] = etabUSER; |
1167 | break; |
1168 | case evdwCUT: |
1169 | tabsel[etiLJ6] = etabLJ6; |
1170 | tabsel[etiLJ12] = etabLJ12; |
1171 | break; |
1172 | case evdwENCADSHIFT: |
1173 | tabsel[etiLJ6] = etabLJ6Encad; |
1174 | tabsel[etiLJ12] = etabLJ12Encad; |
1175 | break; |
1176 | case evdwPME: |
1177 | tabsel[etiLJ6] = etabLJ6Ewald; |
1178 | tabsel[etiLJ12] = etabLJ12; |
1179 | break; |
1180 | default: |
1181 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1181, "Invalid vdwtype %d in %s line %d", vdwtype, |
1182 | __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", __LINE__1182); |
1183 | } |
1184 | |
1185 | if (!b14only && fr->vdw_modifier != eintmodNONE) |
1186 | { |
1187 | if (fr->vdw_modifier != eintmodPOTSHIFT && |
1188 | fr->vdwtype != evdwCUT) |
1189 | { |
1190 | gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off")_gmx_error("incons", "Potential modifiers other than potential-shift are only implemented for LJ cut-off" , "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1190 ); |
1191 | } |
1192 | |
1193 | switch (fr->vdw_modifier) |
1194 | { |
1195 | case eintmodNONE: |
1196 | case eintmodPOTSHIFT: |
1197 | case eintmodEXACTCUTOFF: |
1198 | /* No modification */ |
1199 | break; |
1200 | case eintmodPOTSWITCH: |
1201 | tabsel[etiLJ6] = etabLJ6Switch; |
1202 | tabsel[etiLJ12] = etabLJ12Switch; |
1203 | break; |
1204 | case eintmodFORCESWITCH: |
1205 | tabsel[etiLJ6] = etabLJ6Shift; |
1206 | tabsel[etiLJ12] = etabLJ12Shift; |
1207 | break; |
1208 | default: |
1209 | gmx_incons("Unsupported vdw_modifier")_gmx_error("incons", "Unsupported vdw_modifier", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1209); |
1210 | } |
1211 | } |
1212 | } |
1213 | } |
1214 | |
1215 | t_forcetable make_tables(FILE *out, const output_env_t oenv, |
1216 | const t_forcerec *fr, |
1217 | gmx_bool bVerbose, const char *fn, |
1218 | real rtab, int flags) |
1219 | { |
1220 | const char *fns[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" }; |
1221 | const char *fns14[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" }; |
1222 | FILE *fp; |
1223 | t_tabledata *td; |
1224 | gmx_bool b14only, bReadTab, bGenTab; |
1225 | real x0, y0, yp; |
1226 | int i, j, k, nx, nx0, tabsel[etiNR]; |
1227 | real scalefactor; |
1228 | |
1229 | t_forcetable table; |
1230 | |
1231 | b14only = (flags & GMX_MAKETABLES_14ONLY(1<<1)); |
1232 | |
1233 | if (flags & GMX_MAKETABLES_FORCEUSER(1<<0)) |
1234 | { |
1235 | tabsel[etiCOUL] = etabUSER; |
1236 | tabsel[etiLJ6] = etabUSER; |
1237 | tabsel[etiLJ12] = etabUSER; |
1238 | } |
1239 | else |
1240 | { |
1241 | set_table_type(tabsel, fr, b14only); |
1242 | } |
1243 | snew(td, etiNR)(td) = save_calloc("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1243, (etiNR), sizeof(*(td))); |
1244 | table.r = rtab; |
1245 | table.scale = 0; |
1246 | table.n = 0; |
1247 | table.scale_exp = 0; |
1248 | nx0 = 10; |
1249 | nx = 0; |
1250 | |
1251 | table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP; |
1252 | table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH; |
1253 | table.formatsize = 4; |
1254 | table.ninteractions = 3; |
1255 | table.stride = table.formatsize*table.ninteractions; |
1256 | |
1257 | /* Check whether we have to read or generate */ |
1258 | bReadTab = FALSE0; |
1259 | bGenTab = FALSE0; |
1260 | for (i = 0; (i < etiNR); i++) |
1261 | { |
1262 | if (ETAB_USER(tabsel[i])((tabsel[i]) == etabUSER || (tabsel[i]) == etabEwaldUser || ( tabsel[i]) == etabEwaldUserSwitch)) |
1263 | { |
1264 | bReadTab = TRUE1; |
1265 | } |
1266 | if (tabsel[i] != etabUSER) |
1267 | { |
1268 | bGenTab = TRUE1; |
1269 | } |
1270 | } |
1271 | if (bReadTab) |
1272 | { |
1273 | read_tables(out, fn, etiNR, 0, td); |
1274 | if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY(1<<1))) |
1275 | { |
1276 | rtab = td[0].x[td[0].nx-1]; |
1277 | table.n = td[0].nx; |
1278 | nx = table.n; |
1279 | } |
1280 | else |
1281 | { |
1282 | if (td[0].x[td[0].nx-1] < rtab) |
1283 | { |
1284 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1284, "Tables in file %s not long enough for cut-off:\n" |
1285 | "\tshould be at least %f nm\n", fn, rtab); |
1286 | } |
1287 | nx = table.n = (int)(rtab*td[0].tabscale + 0.5); |
1288 | } |
1289 | table.scale = td[0].tabscale; |
1290 | nx0 = td[0].nx0; |
1291 | } |
1292 | if (bGenTab) |
1293 | { |
1294 | if (!bReadTab) |
1295 | { |
1296 | #ifdef GMX_DOUBLE |
1297 | table.scale = 2000.0; |
1298 | #else |
1299 | table.scale = 500.0; |
1300 | #endif |
1301 | nx = table.n = rtab*table.scale; |
1302 | } |
1303 | } |
1304 | if (fr->bBHAM) |
1305 | { |
1306 | if (fr->bham_b_max != 0) |
1307 | { |
1308 | table.scale_exp = table.scale/fr->bham_b_max; |
1309 | } |
1310 | else |
1311 | { |
1312 | table.scale_exp = table.scale; |
1313 | } |
1314 | } |
1315 | |
1316 | /* Each table type (e.g. coul,lj6,lj12) requires four |
1317 | * numbers per nx+1 data points. For performance reasons we want |
1318 | * the table data to be aligned to 16-byte. |
1319 | */ |
1320 | snew_aligned(table.data, 12*(nx+1)*sizeof(real), 32)(table.data) = save_calloc_aligned("table.data", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1320, (12*(nx+1)*sizeof(real)), sizeof(*(table.data)), 32); |
1321 | |
1322 | for (k = 0; (k < etiNR); k++) |
1323 | { |
1324 | if (tabsel[k] != etabUSER) |
1325 | { |
1326 | init_table(nx, nx0, |
1327 | (tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale, |
1328 | &(td[k]), !bReadTab); |
1329 | fill_table(&(td[k]), tabsel[k], fr); |
1330 | if (out) |
1331 | { |
1332 | fprintf(out, "%s table with %d data points for %s%s.\n" |
1333 | "Tabscale = %g points/nm\n", |
1334 | ETAB_USER(tabsel[k])((tabsel[k]) == etabUSER || (tabsel[k]) == etabEwaldUser || ( tabsel[k]) == etabEwaldUserSwitch) ? "Modified" : "Generated", |
1335 | td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name, |
1336 | td[k].tabscale); |
1337 | } |
1338 | } |
1339 | |
1340 | /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels |
1341 | * by including the derivative constants (6.0 or 12.0) in the parameters, since |
1342 | * we no longer calculate force in most steps. This means the c6/c12 parameters |
1343 | * have been scaled up, so we need to scale down the table interactions too. |
1344 | * It comes here since we need to scale user tables too. |
1345 | */ |
1346 | if (k == etiLJ6) |
1347 | { |
1348 | scalefactor = 1.0/6.0; |
1349 | } |
1350 | else if (k == etiLJ12 && tabsel[k] != etabEXPMIN) |
1351 | { |
1352 | scalefactor = 1.0/12.0; |
1353 | } |
1354 | else |
1355 | { |
1356 | scalefactor = 1.0; |
1357 | } |
1358 | |
1359 | copy2table(table.n, k*4, 12, td[k].x, td[k].v, td[k].f, scalefactor, table.data); |
1360 | |
1361 | if (bDebugMode() && bVerbose) |
1362 | { |
1363 | if (b14only) |
1364 | { |
1365 | fp = xvgropen(fns14[k], fns14[k], "r", "V", oenv); |
1366 | } |
1367 | else |
1368 | { |
1369 | fp = xvgropen(fns[k], fns[k], "r", "V", oenv); |
1370 | } |
1371 | /* plot the output 5 times denser than the table data */ |
1372 | for (i = 5*((nx0+1)/2); i < 5*table.n; i++) |
1373 | { |
1374 | x0 = i*table.r/(5*(table.n-1)); |
1375 | evaluate_table(table.data, 4*k, 12, table.scale, x0, &y0, &yp); |
1376 | fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp); |
1377 | } |
1378 | gmx_fio_fclose(fp); |
1379 | } |
1380 | done_tabledata(&(td[k])); |
1381 | } |
1382 | sfree(td)save_free("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1382, (td)); |
1383 | |
1384 | return table; |
1385 | } |
1386 | |
1387 | t_forcetable make_gb_table(const output_env_t oenv, |
1388 | const t_forcerec *fr) |
1389 | { |
1390 | const char *fns[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" }; |
1391 | const char *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" }; |
1392 | FILE *fp; |
1393 | t_tabledata *td; |
1394 | gmx_bool bReadTab, bGenTab; |
1395 | real x0, y0, yp; |
1396 | int i, j, k, nx, nx0, tabsel[etiNR]; |
1397 | double r, r2, Vtab, Ftab, expterm; |
1398 | |
1399 | t_forcetable table; |
1400 | |
1401 | double abs_error_r, abs_error_r2; |
1402 | double rel_error_r, rel_error_r2; |
1403 | double rel_error_r_old = 0, rel_error_r2_old = 0; |
1404 | double x0_r_error, x0_r2_error; |
1405 | |
1406 | |
1407 | /* Only set a Coulomb table for GB */ |
1408 | /* |
1409 | tabsel[0]=etabGB; |
1410 | tabsel[1]=-1; |
1411 | tabsel[2]=-1; |
1412 | */ |
1413 | |
1414 | /* Set the table dimensions for GB, not really necessary to |
1415 | * use etiNR (since we only have one table, but ...) |
1416 | */ |
1417 | snew(td, 1)(td) = save_calloc("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1417, (1), sizeof(*(td))); |
1418 | table.interaction = GMX_TABLE_INTERACTION_ELEC; |
1419 | table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH; |
1420 | table.r = fr->gbtabr; |
1421 | table.scale = fr->gbtabscale; |
1422 | table.scale_exp = 0; |
1423 | table.n = table.scale*table.r; |
1424 | table.formatsize = 4; |
1425 | table.ninteractions = 1; |
1426 | table.stride = table.formatsize*table.ninteractions; |
1427 | nx0 = 0; |
1428 | nx = table.scale*table.r; |
1429 | |
1430 | /* Check whether we have to read or generate |
1431 | * We will always generate a table, so remove the read code |
1432 | * (Compare with original make_table function |
1433 | */ |
1434 | bReadTab = FALSE0; |
1435 | bGenTab = TRUE1; |
Value stored to 'bGenTab' is never read | |
1436 | |
1437 | /* Each table type (e.g. coul,lj6,lj12) requires four |
1438 | * numbers per datapoint. For performance reasons we want |
1439 | * the table data to be aligned to 16-byte. This is accomplished |
1440 | * by allocating 16 bytes extra to a temporary pointer, and then |
1441 | * calculating an aligned pointer. This new pointer must not be |
1442 | * used in a free() call, but thankfully we're sloppy enough not |
1443 | * to do this :-) |
1444 | */ |
1445 | |
1446 | snew_aligned(table.data, 4*nx, 32)(table.data) = save_calloc_aligned("table.data", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1446, (4*nx), sizeof(*(table.data)), 32); |
1447 | |
1448 | init_table(nx, nx0, table.scale, &(td[0]), !bReadTab); |
1449 | |
1450 | /* Local implementation so we don't have to use the etabGB |
1451 | * enum above, which will cause problems later when |
1452 | * making the other tables (right now even though we are using |
1453 | * GB, the normal Coulomb tables will be created, but this |
1454 | * will cause a problem since fr->eeltype==etabGB which will not |
1455 | * be defined in fill_table and set_table_type |
1456 | */ |
1457 | |
1458 | for (i = nx0; i < nx; i++) |
1459 | { |
1460 | r = td->x[i]; |
1461 | r2 = r*r; |
1462 | expterm = exp(-0.25*r2); |
1463 | |
1464 | Vtab = 1/sqrt(r2+expterm); |
1465 | Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm)); |
1466 | |
1467 | /* Convert to single precision when we store to mem */ |
1468 | td->v[i] = Vtab; |
1469 | td->f[i] = Ftab; |
1470 | |
1471 | } |
1472 | |
1473 | copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data); |
1474 | |
1475 | if (bDebugMode()) |
1476 | { |
1477 | fp = xvgropen(fns[0], fns[0], "r", "V", oenv); |
1478 | /* plot the output 5 times denser than the table data */ |
1479 | /* for(i=5*nx0;i<5*table.n;i++) */ |
1480 | for (i = nx0; i < table.n; i++) |
1481 | { |
1482 | /* x0=i*table.r/(5*table.n); */ |
1483 | x0 = i*table.r/table.n; |
1484 | evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp); |
1485 | fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp); |
1486 | |
1487 | } |
1488 | gmx_fio_fclose(fp); |
1489 | } |
1490 | |
1491 | /* |
1492 | for(i=100*nx0;i<99.81*table.n;i++) |
1493 | { |
1494 | r = i*table.r/(100*table.n); |
1495 | r2 = r*r; |
1496 | expterm = exp(-0.25*r2); |
1497 | |
1498 | Vtab = 1/sqrt(r2+expterm); |
1499 | Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm)); |
1500 | |
1501 | |
1502 | evaluate_table(table.data,0,4,table.scale,r,&y0,&yp); |
1503 | printf("gb: i=%d, x0=%g, y0=%15.15f, Vtab=%15.15f, yp=%15.15f, Ftab=%15.15f\n",i,r, y0, Vtab, yp, Ftab); |
1504 | |
1505 | abs_error_r=fabs(y0-Vtab); |
1506 | abs_error_r2=fabs(yp-(-1)*Ftab); |
1507 | |
1508 | rel_error_r=abs_error_r/y0; |
1509 | rel_error_r2=fabs(abs_error_r2/yp); |
1510 | |
1511 | |
1512 | if(rel_error_r>rel_error_r_old) |
1513 | { |
1514 | rel_error_r_old=rel_error_r; |
1515 | x0_r_error=x0; |
1516 | } |
1517 | |
1518 | if(rel_error_r2>rel_error_r2_old) |
1519 | { |
1520 | rel_error_r2_old=rel_error_r2; |
1521 | x0_r2_error=x0; |
1522 | } |
1523 | } |
1524 | |
1525 | printf("gb: MAX REL ERROR IN R=%15.15f, MAX REL ERROR IN R2=%15.15f\n",rel_error_r_old, rel_error_r2_old); |
1526 | printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error); |
1527 | |
1528 | exit(1); */ |
1529 | done_tabledata(&(td[0])); |
1530 | sfree(td)save_free("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1530, (td)); |
1531 | |
1532 | return table; |
1533 | |
1534 | |
1535 | } |
1536 | |
1537 | t_forcetable make_atf_table(FILE *out, const output_env_t oenv, |
1538 | const t_forcerec *fr, |
1539 | const char *fn, |
1540 | matrix box) |
1541 | { |
1542 | const char *fns[3] = { "tf_tab.xvg", "atfdtab.xvg", "atfrtab.xvg" }; |
1543 | FILE *fp; |
1544 | t_tabledata *td; |
1545 | real x0, y0, yp, rtab; |
1546 | int i, nx, nx0; |
1547 | real rx, ry, rz, box_r; |
1548 | |
1549 | t_forcetable table; |
1550 | |
1551 | |
1552 | /* Set the table dimensions for ATF, not really necessary to |
1553 | * use etiNR (since we only have one table, but ...) |
1554 | */ |
1555 | snew(td, 1)(td) = save_calloc("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1555, (1), sizeof(*(td))); |
1556 | |
1557 | if (fr->adress_type == eAdressSphere) |
1558 | { |
1559 | /* take half box diagonal direction as tab range */ |
1560 | rx = 0.5*box[0][0]+0.5*box[1][0]+0.5*box[2][0]; |
1561 | ry = 0.5*box[0][1]+0.5*box[1][1]+0.5*box[2][1]; |
1562 | rz = 0.5*box[0][2]+0.5*box[1][2]+0.5*box[2][2]; |
1563 | box_r = sqrt(rx*rx+ry*ry+rz*rz); |
1564 | |
1565 | } |
1566 | else |
1567 | { |
1568 | /* xsplit: take half box x direction as tab range */ |
1569 | box_r = box[0][0]/2; |
1570 | } |
1571 | table.r = box_r; |
1572 | table.scale = 0; |
1573 | table.n = 0; |
1574 | table.scale_exp = 0; |
1575 | nx0 = 10; |
1576 | nx = 0; |
1577 | |
1578 | read_tables(out, fn, 1, 0, td); |
1579 | rtab = td[0].x[td[0].nx-1]; |
1580 | |
1581 | if (fr->adress_type == eAdressXSplit && (rtab < box[0][0]/2)) |
1582 | { |
1583 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1583, "AdResS full box therm force table in file %s extends to %f:\n" |
1584 | "\tshould extend to at least half the length of the box in x-direction" |
1585 | "%f\n", fn, rtab, box[0][0]/2); |
1586 | } |
1587 | if (rtab < box_r) |
1588 | { |
1589 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c", 1589, "AdResS full box therm force table in file %s extends to %f:\n" |
1590 | "\tshould extend to at least for spherical adress" |
1591 | "%f (=distance from center to furthermost point in box \n", fn, rtab, box_r); |
1592 | } |
1593 | |
1594 | |
1595 | table.n = td[0].nx; |
1596 | nx = table.n; |
1597 | table.scale = td[0].tabscale; |
1598 | nx0 = td[0].nx0; |
1599 | |
1600 | /* Each table type (e.g. coul,lj6,lj12) requires four |
1601 | * numbers per datapoint. For performance reasons we want |
1602 | * the table data to be aligned to 16-byte. This is accomplished |
1603 | * by allocating 16 bytes extra to a temporary pointer, and then |
1604 | * calculating an aligned pointer. This new pointer must not be |
1605 | * used in a free() call, but thankfully we're sloppy enough not |
1606 | * to do this :-) |
1607 | */ |
1608 | |
1609 | snew_aligned(table.data, 4*nx, 32)(table.data) = save_calloc_aligned("table.data", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1609, (4*nx), sizeof(*(table.data)), 32); |
1610 | |
1611 | copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data); |
1612 | |
1613 | if (bDebugMode()) |
1614 | { |
1615 | fp = xvgropen(fns[0], fns[0], "r", "V", oenv); |
1616 | /* plot the output 5 times denser than the table data */ |
1617 | /* for(i=5*nx0;i<5*table.n;i++) */ |
1618 | |
1619 | for (i = 5*((nx0+1)/2); i < 5*table.n; i++) |
1620 | { |
1621 | /* x0=i*table.r/(5*table.n); */ |
1622 | x0 = i*table.r/(5*(table.n-1)); |
1623 | evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp); |
1624 | fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp); |
1625 | |
1626 | } |
1627 | gmx_ffclose(fp); |
1628 | } |
1629 | |
1630 | done_tabledata(&(td[0])); |
1631 | sfree(td)save_free("td", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1631, (td)); |
1632 | |
1633 | table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP; |
1634 | table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH; |
1635 | table.formatsize = 4; |
1636 | table.ninteractions = 3; |
1637 | table.stride = table.formatsize*table.ninteractions; |
1638 | |
1639 | |
1640 | return table; |
1641 | } |
1642 | |
1643 | bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle) |
1644 | { |
1645 | t_tabledata td; |
1646 | double start; |
1647 | int i; |
1648 | bondedtable_t tab; |
1649 | |
1650 | if (angle < 2) |
1651 | { |
1652 | start = 0; |
1653 | } |
1654 | else |
1655 | { |
1656 | start = -180.0; |
1657 | } |
1658 | read_tables(fplog, fn, 1, angle, &td); |
1659 | if (angle > 0) |
1660 | { |
1661 | /* Convert the table from degrees to radians */ |
1662 | for (i = 0; i < td.nx; i++) |
1663 | { |
1664 | td.x[i] *= DEG2RAD(3.14159265358979323846/180.0); |
1665 | td.f[i] *= RAD2DEG(180.0/3.14159265358979323846); |
1666 | } |
1667 | td.tabscale *= RAD2DEG(180.0/3.14159265358979323846); |
1668 | } |
1669 | tab.n = td.nx; |
1670 | tab.scale = td.tabscale; |
1671 | snew(tab.data, tab.n*4)(tab.data) = save_calloc("tab.data", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/tables.c" , 1671, (tab.n*4), sizeof(*(tab.data))); |
1672 | copy2table(tab.n, 0, 4, td.x, td.v, td.f, 1.0, tab.data); |
1673 | done_tabledata(&td); |
1674 | |
1675 | return tab; |
1676 | } |