Bug Summary

File:gromacs/mdlib/tables.c
Location:line 1391, column 21
Description:Value stored to 'fns14' during its initialization is never read

Annotated Source Code

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 */
58enum {
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
87typedef 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 */
94static 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 */
119enum {
120 etiCOUL, etiLJ6, etiLJ12, etiNR
121};
122
123typedef 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
134double 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
146double 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
164void 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 */
316real 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 */
344static 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
365static 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
403static 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
423static 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
523static 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
563static 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
734static 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
748static 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
1054static 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
1215t_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
1387t_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" };
Value stored to 'fns14' during its initialization is never read
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;
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
1537t_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
1643bondedtable_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}