3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
44 #include "gmx_fatal.h"
54 /* All the possible (implemented) table functions */
79 /** Evaluates to true if the table type contains user data. */
80 #define ETAB_USER(e) ((e) == etabUSER || \
81 (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
88 /* This structure holds name and a flag that tells whether
89 this is a Coulomb type funtion */
90 static const t_tab_props tprops[etabNR] = {
93 { "LJ6Shift", FALSE },
94 { "LJ12Shift", FALSE },
100 { "Ewald-Switch", TRUE },
101 { "Ewald-User", TRUE },
102 { "Ewald-User-Switch", TRUE },
103 { "LJ6Switch", FALSE },
104 { "LJ12Switch", FALSE },
105 { "COULSwitch", TRUE },
106 { "LJ6-Encad shift", FALSE },
107 { "LJ12-Encad shift", FALSE },
108 { "COUL-Encad shift", TRUE },
113 /* Index in the table that says which function to use */
114 enum { etiCOUL, etiLJ6, etiLJ12, etiNR };
122 #define pow2(x) ((x)*(x))
123 #define pow3(x) ((x)*(x)*(x))
124 #define pow4(x) ((x)*(x)*(x)*(x))
125 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
127 /* Calculate the potential and force for an r value
128 * in exactly the same way it is done in the inner loop.
129 * VFtab is a pointer to the table data, offset is
130 * the point where we should begin and stride is
131 * 4 if we have a buckingham table, 3 otherwise.
132 * If you want to evaluate table no N, set offset to 4*N.
134 * We use normal precision here, since that is what we
135 * will use in the inner loops.
137 static void evaluate_table(real VFtab[], int offset, int stride,
138 real tabscale, real r, real *y, real *yp)
142 real Y,F,Geps,Heps2,Fp;
151 Geps = eps*VFtab[n+2];
152 Heps2 = eps2*VFtab[n+3];
155 *yp = (Fp+Geps+2.0*Heps2)*tabscale;
159 static void splint(real xa[],real ya[],real y2a[],
160 int n,real x,real *y,real *yp)
169 while ((khi-klo) > 1) {
178 gmx_fatal(FARGS,"Bad XA input to function splint");
181 *y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
182 *yp = (ya[khi]-ya[klo])/h+((3*a*a-1)*y2a[klo]-(3*b*b-1)*y2a[khi])*h/6.0;
185 F = (ya[khi]-ya[klo]-(h*h/6.0)*(2*y2a[klo]+y2a[khi]));
186 G = (h*h/2.0)*y2a[klo];
187 H = (h*h/6.0)*(y2a[khi]-y2a[klo]);
188 *y = ya[klo] + eps*F + eps*eps*G + eps*eps*eps*H;
189 *yp = (F + 2*eps*G + 3*eps*eps*H)/h;
193 static void copy2table(int n,int offset,int stride,
194 double x[],double Vtab[],double Ftab[],
197 /* Use double prec. for the intermediary variables
198 * and temporary x/vtab/vtab2 data to avoid unnecessary
205 for(i=0; (i<n); i++) {
209 G = 3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
210 H = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] + Ftab[i])*h;
212 /* Fill the last entry with a linear potential,
213 * this is mainly for rounding issues with angle and dihedral potentials.
219 nn0 = offset + i*stride;
227 static void init_table(FILE *fp,int n,int nx0,
228 double tabscale,t_tabledata *td,gmx_bool bAlloc)
234 td->tabscale = tabscale;
240 for(i=0; (i<td->nx); i++)
241 td->x[i] = i/tabscale;
244 static void spline_forces(int nx,double h,double v[],gmx_bool bS3,gmx_bool bE3,
251 /* Formulas can be found in:
252 * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
255 if (nx < 4 && (bS3 || bE3))
256 gmx_fatal(FARGS,"Can not generate splines with third derivative boundary conditions with less than 4 (%d) points",nx);
258 /* To make life easy we initially set the spacing to 1
259 * and correct for this at the end.
263 /* Fit V''' at the start */
264 v3 = v[3] - 3*v[2] + 3*v[1] - v[0];
266 fprintf(debug,"The left third derivative is %g\n",v3/(h*h*h));
267 b_s = 2*(v[1] - v[0]) + v3/6;
271 /* Fit V'' at the start */
274 v2 = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
275 /* v2 = v[2] - 2*v[1] + v[0]; */
277 fprintf(debug,"The left second derivative is %g\n",v2/(h*h));
278 b_s = 3*(v[1] - v[0]) - v2/2;
282 b_s = 3*(v[2] - v[0]) + f[0]*h;
286 /* Fit V''' at the end */
287 v3 = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
289 fprintf(debug,"The right third derivative is %g\n",v3/(h*h*h));
290 b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
293 /* V'=0 at the end */
294 b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
299 beta = (bS3 ? 1 : 4);
301 /* For V'' fitting */
302 /* beta = (bS3 ? 2 : 4); */
305 for(i=start+1; i<end; i++) {
308 b = 3*(v[i+1] - v[i-1]);
309 f[i] = (b - f[i-1])/beta;
311 gamma[end-1] = 1/beta;
312 beta = (bE3 ? 1 : 4) - gamma[end-1];
313 f[end-1] = (b_e - f[end-2])/beta;
315 for(i=end-2; i>=start; i--)
316 f[i] -= gamma[i+1]*f[i+1];
319 /* Correct for the minus sign and the spacing */
320 for(i=start; i<end; i++)
324 static void set_forces(FILE *fp,int angle,
325 int nx,double h,double v[],double f[],
332 "Force generation for dihedral tables is not (yet) implemented");
335 while (v[start] == 0)
347 fprintf(fp,"Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
348 table+1,start*h,end==nx ? "V'''" : "V'=0",(end-1)*h);
349 spline_forces(end-start,h,v+start,TRUE,end==nx,f+start);
352 static void read_tables(FILE *fp,const char *fn,
353 int ntab,int angle,t_tabledata td[])
357 double **yy=NULL,start,end,dx0,dx1,ssd,vm,vp,f,numf;
358 int k,i,nx,nx0=0,ny,nny,ns;
359 gmx_bool bAllZero,bZeroV,bZeroF;
363 libfn = gmxlibfn(fn);
364 nx = read_xvg(libfn,&yy,&ny);
366 gmx_fatal(FARGS,"Trying to read file %s, but nr columns = %d, should be %d",
371 "The first distance in file %s is %f nm instead of %f nm",
379 if (yy[0][0] != start || yy[0][nx-1] != end)
380 gmx_fatal(FARGS,"The angles in file %s should go from %f to %f instead of %f to %f\n",
381 libfn,start,end,yy[0][0],yy[0][nx-1]);
384 tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
387 fprintf(fp,"Read user tables from %s with %d data points.\n",libfn,nx);
389 fprintf(fp,"Tabscale = %g points/nm\n",tabscale);
393 for(k=0; k<ntab; k++) {
396 for(i=0; (i < nx); i++) {
398 dx0 = yy[0][i-1] - yy[0][i-2];
399 dx1 = yy[0][i] - yy[0][i-1];
400 /* Check for 1% deviation in spacing */
401 if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1))) {
402 gmx_fatal(FARGS,"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]);
405 if (yy[1+k*2][i] != 0) {
411 if (yy[1+k*2][i] > 0.01*GMX_REAL_MAX ||
412 yy[1+k*2][i] < -0.01*GMX_REAL_MAX) {
413 gmx_fatal(FARGS,"Out of range potential value %g in file '%s'",
417 if (yy[1+k*2+1][i] != 0) {
423 if (yy[1+k*2+1][i] > 0.01*GMX_REAL_MAX ||
424 yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX) {
425 gmx_fatal(FARGS,"Out of range force value %g in file '%s'",
431 if (!bZeroV && bZeroF) {
432 set_forces(fp,angle,nx,1/tabscale,yy[1+k*2],yy[1+k*2+1],k);
434 /* Check if the second column is close to minus the numerical
435 * derivative of the first column.
439 for(i=1; (i < nx-1); i++) {
443 if (vm != 0 && vp != 0 && f != 0) {
444 /* Take the centered difference */
445 numf = -(vp - vm)*0.5*tabscale;
446 ssd += fabs(2*(f - numf)/(f + numf));
452 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));
454 fprintf(debug,"%s",buf);
457 fprintf(fp,"\nWARNING: %s\n",buf);
458 fprintf(stderr,"\nWARNING: %s\n",buf);
463 if (bAllZero && fp) {
464 fprintf(fp,"\nNOTE: All elements in table %s are zero\n\n",libfn);
467 for(k=0; (k<ntab); k++) {
468 init_table(fp,nx,nx0,tabscale,&(td[k]),TRUE);
469 for(i=0; (i<nx); i++) {
470 td[k].x[i] = yy[0][i];
471 td[k].v[i] = yy[2*k+1][i];
472 td[k].f[i] = yy[2*k+2][i];
475 for(i=0; (i<ny); i++)
481 static void done_tabledata(t_tabledata *td)
493 static void fill_table(t_tabledata *td,int tp,const t_forcerec *fr)
495 /* Fill the table according to the formulas in the manual.
496 * In principle, we only need the potential and the second
497 * derivative, but then we would have to do lots of calculations
498 * in the inner loop. By precalculating some terms (see manual)
499 * we get better eventual performance, despite a larger table.
501 * Since some of these higher-order terms are very small,
502 * we always use double precision to calculate them here, in order
503 * to avoid unnecessary loss of precision.
510 double r1,rc,r12,r13;
512 double expr,Vtab,Ftab;
513 /* Parameters for David's function */
514 double A=0,B=0,C=0,A_3=0,B_4=0;
515 /* Parameters for the switching function */
517 /* Temporary parameters */
518 gmx_bool bSwitch,bShift;
519 double ewc=fr->ewaldcoeff;
520 double isp= 0.564189583547756;
522 bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
523 (tp == etabCOULSwitch) ||
524 (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
525 bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
530 if (tprops[tp].bCoulomb) {
531 r1 = fr->rcoulomb_switch;
535 r1 = fr->rvdw_switch;
539 ksw = 1.0/(pow5(rc-r1));
545 else if (tp == etabLJ6Shift)
550 A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc,p+2)*pow2(rc-r1));
551 B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc,p+2)*pow3(rc-r1));
552 C = 1.0/pow(rc,p)-A/3.0*pow3(rc-r1)-B/4.0*pow4(rc-r1);
553 if (tp == etabLJ6Shift) {
561 if (debug) { fprintf(debug,"Setting up tables\n"); fflush(debug); }
564 fp=xvgropen("switch.xvg","switch","r","s");
567 for(i=td->nx0; (i<td->nx); i++) {
571 if (gmx_within_tol(reppow,12.0,10*GMX_DOUBLE_EPS)) {
574 r12 = pow(r,-reppow);
579 /* swi is function, swi1 1st derivative and swi2 2nd derivative */
580 /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
581 * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
583 * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
592 swi = 1 - 10*pow3(r-r1)*ksw*pow2(rc-r1)
593 + 15*pow4(r-r1)*ksw*(rc-r1) - 6*pow5(r-r1)*ksw;
594 swi1 = -30*pow2(r-r1)*ksw*pow2(rc-r1)
595 + 60*pow3(r-r1)*ksw*(rc-r1) - 30*pow4(r-r1)*ksw;
598 else { /* not really needed, but avoids compiler warnings... */
603 fprintf(fp,"%10g %10g %10g %10g\n",r,swi,swi1,swi2);
626 Ftab = reppow*Vtab/r;
633 Ftab = reppow*Vtab/r;
638 Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
639 Ftab = -(6.0*r6/r-6.0*rc6/rc);
647 Vtab = r12-12.0*(rc-r)*rc6*rc6/rc-1.0*rc6*rc6;
648 Ftab = 12.0*r12/r-12.0*rc6*rc6/rc;
666 case etabEwaldSwitch:
667 Vtab = gmx_erfc(ewc*r)/r;
668 Ftab = gmx_erfc(ewc*r)/r2+2*exp(-(ewc*ewc*r2))*ewc*isp/r;
671 case etabEwaldUserSwitch:
672 /* Only calculate minus the reciprocal space contribution */
673 Vtab = -gmx_erf(ewc*r)/r;
674 Ftab = -gmx_erf(ewc*r)/r2+2*exp(-(ewc*ewc*r2))*ewc*isp/r;
678 Vtab = 1.0/r + fr->k_rf*r2 - fr->c_rf;
679 Ftab = 1.0/r2 - 2*fr->k_rf*r;
680 if (tp == etabRF_ZERO && r >= rc) {
692 Vtab = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
693 Ftab = 1.0/r2-1.0/(rc*rc);
700 gmx_fatal(FARGS,"Table type %d not implemented yet. (%s,%d)",
701 tp,__FILE__,__LINE__);
704 /* Normal coulomb with cut-off correction for potential */
707 /* If in Shifting range add something to it */
711 Vtab += - A_3*r13 - B_4*r12*r12;
712 Ftab += A*r12 + B*r13;
722 if ((r > r1) && bSwitch) {
723 Ftab = Ftab*swi - Vtab*swi1;
727 /* Convert to single precision when we store to mem */
732 /* Continue the table linearly from nx0 to 0.
733 * These values are only required for energy minimization with overlap or TPI.
735 for(i=td->nx0-1; i>=0; i--) {
736 td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
737 td->f[i] = td->f[i+1];
745 static void set_table_type(int tabsel[],const t_forcerec *fr,gmx_bool b14only)
749 /* Set the different table indices.
755 switch (fr->eeltype) {
761 case eelPMEUSERSWITCH:
768 eltype = fr->eeltype;
773 tabsel[etiCOUL] = etabCOUL;
777 tabsel[etiCOUL] = etabShift;
780 if (fr->rcoulomb > fr->rcoulomb_switch)
781 tabsel[etiCOUL] = etabShift;
783 tabsel[etiCOUL] = etabCOUL;
787 tabsel[etiCOUL] = etabEwald;
790 tabsel[etiCOUL] = etabEwaldSwitch;
793 tabsel[etiCOUL] = etabEwaldUser;
795 case eelPMEUSERSWITCH:
796 tabsel[etiCOUL] = etabEwaldUserSwitch;
801 tabsel[etiCOUL] = etabRF;
804 tabsel[etiCOUL] = etabRF_ZERO;
807 tabsel[etiCOUL] = etabCOULSwitch;
810 tabsel[etiCOUL] = etabUSER;
813 tabsel[etiCOUL] = etabCOULEncad;
816 gmx_fatal(FARGS,"Invalid eeltype %d",eltype);
819 /* Van der Waals time */
820 if (fr->bBHAM && !b14only) {
821 tabsel[etiLJ6] = etabLJ6;
822 tabsel[etiLJ12] = etabEXPMIN;
824 if (b14only && fr->vdwtype != evdwUSER)
827 vdwtype = fr->vdwtype;
831 tabsel[etiLJ6] = etabLJ6Switch;
832 tabsel[etiLJ12] = etabLJ12Switch;
835 tabsel[etiLJ6] = etabLJ6Shift;
836 tabsel[etiLJ12] = etabLJ12Shift;
839 tabsel[etiLJ6] = etabUSER;
840 tabsel[etiLJ12] = etabUSER;
843 tabsel[etiLJ6] = etabLJ6;
844 tabsel[etiLJ12] = etabLJ12;
847 tabsel[etiLJ6] = etabLJ6Encad;
848 tabsel[etiLJ12] = etabLJ12Encad;
851 gmx_fatal(FARGS,"Invalid vdwtype %d in %s line %d",vdwtype,
857 t_forcetable make_tables(FILE *out,const output_env_t oenv,
858 const t_forcerec *fr,
859 gmx_bool bVerbose,const char *fn,
862 const char *fns[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
863 const char *fns14[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" };
866 gmx_bool b14only,bReadTab,bGenTab;
868 int i,j,k,nx,nx0,tabsel[etiNR];
872 b14only = (flags & GMX_MAKETABLES_14ONLY);
874 if (flags & GMX_MAKETABLES_FORCEUSER) {
875 tabsel[etiCOUL] = etabUSER;
876 tabsel[etiLJ6] = etabUSER;
877 tabsel[etiLJ12] = etabUSER;
879 set_table_type(tabsel,fr,b14only);
889 /* Check whether we have to read or generate */
892 for(i=0; (i<etiNR); i++) {
893 if (ETAB_USER(tabsel[i]))
895 if (tabsel[i] != etabUSER)
899 read_tables(out,fn,etiNR,0,td);
900 if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY)) {
901 rtab = td[0].x[td[0].nx-1];
905 if (td[0].x[td[0].nx-1] < rtab)
906 gmx_fatal(FARGS,"Tables in file %s not long enough for cut-off:\n"
907 "\tshould be at least %f nm\n",fn,rtab);
908 nx = table.n = (int)(rtab*td[0].tabscale + 0.5);
910 table.scale = td[0].tabscale;
916 table.scale = 2000.0;
920 nx = table.n = rtab*table.scale;
924 if(fr->bham_b_max!=0)
925 table.scale_exp = table.scale/fr->bham_b_max;
927 table.scale_exp = table.scale;
930 /* Each table type (e.g. coul,lj6,lj12) requires four
931 * numbers per nx+1 data points. For performance reasons we want
932 * the table data to be aligned to 16-byte.
934 snew_aligned(table.tab, 12*(nx+1)*sizeof(real),16);
936 for(k=0; (k<etiNR); k++) {
937 if (tabsel[k] != etabUSER) {
938 init_table(out,nx,nx0,
939 (tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
941 fill_table(&(td[k]),tabsel[k],fr);
943 fprintf(out,"%s table with %d data points for %s%s.\n"
944 "Tabscale = %g points/nm\n",
945 ETAB_USER(tabsel[k]) ? "Modified" : "Generated",
946 td[k].nx,b14only?"1-4 ":"",tprops[tabsel[k]].name,
949 copy2table(table.n,k*4,12,td[k].x,td[k].v,td[k].f,table.tab);
951 if (bDebugMode() && bVerbose) {
953 fp=xvgropen(fns14[k],fns14[k],"r","V",oenv);
955 fp=xvgropen(fns[k],fns[k],"r","V",oenv);
956 /* plot the output 5 times denser than the table data */
957 for(i=5*((nx0+1)/2); i<5*table.n; i++) {
958 x0 = i*table.r/(5*(table.n-1));
959 evaluate_table(table.tab,4*k,12,table.scale,x0,&y0,&yp);
960 fprintf(fp,"%15.10e %15.10e %15.10e\n",x0,y0,yp);
964 done_tabledata(&(td[k]));
971 t_forcetable make_gb_table(FILE *out,const output_env_t oenv,
972 const t_forcerec *fr,
976 const char *fns[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
977 const char *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
980 gmx_bool bReadTab,bGenTab;
982 int i,j,k,nx,nx0,tabsel[etiNR];
984 double r,r2,Vtab,Ftab,expterm;
988 double abs_error_r, abs_error_r2;
989 double rel_error_r, rel_error_r2;
990 double rel_error_r_old=0, rel_error_r2_old=0;
991 double x0_r_error, x0_r2_error;
994 /* Only set a Coulomb table for GB */
1001 /* Set the table dimensions for GB, not really necessary to
1002 * use etiNR (since we only have one table, but ...)
1005 table.r = fr->gbtabr;
1006 table.scale = fr->gbtabscale;
1007 table.scale_exp = 0;
1008 table.n = table.scale*table.r;
1010 nx = table.scale*table.r;
1012 /* Check whether we have to read or generate
1013 * We will always generate a table, so remove the read code
1014 * (Compare with original make_table function
1019 /* Each table type (e.g. coul,lj6,lj12) requires four
1020 * numbers per datapoint. For performance reasons we want
1021 * the table data to be aligned to 16-byte. This is accomplished
1022 * by allocating 16 bytes extra to a temporary pointer, and then
1023 * calculating an aligned pointer. This new pointer must not be
1024 * used in a free() call, but thankfully we're sloppy enough not
1028 /* 4 fp entries per table point, nx+1 points, and 16 bytes extra
1030 p_tmp = malloc(4*(nx+1)*sizeof(real)+16);
1032 /* align it - size_t has the same same as a pointer */
1033 table.tab = (real *) (((size_t) p_tmp + 16) & (~((size_t) 15)));
1035 init_table(out,nx,nx0,table.scale,&(td[0]),!bReadTab);
1037 /* Local implementation so we don't have to use the etabGB
1038 * enum above, which will cause problems later when
1039 * making the other tables (right now even though we are using
1040 * GB, the normal Coulomb tables will be created, but this
1041 * will cause a problem since fr->eeltype==etabGB which will not
1042 * be defined in fill_table and set_table_type
1051 expterm = exp(-0.25*r2);
1053 Vtab = 1/sqrt(r2+expterm);
1054 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1056 /* Convert to single precision when we store to mem */
1062 copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,table.tab);
1066 fp=xvgropen(fns[0],fns[0],"r","V",oenv);
1067 /* plot the output 5 times denser than the table data */
1068 /* for(i=5*nx0;i<5*table.n;i++) */
1069 for(i=nx0;i<table.n;i++)
1071 /* x0=i*table.r/(5*table.n); */
1072 x0=i*table.r/table.n;
1073 evaluate_table(table.tab,0,4,table.scale,x0,&y0,&yp);
1074 fprintf(fp,"%15.10e %15.10e %15.10e\n",x0,y0,yp);
1081 for(i=100*nx0;i<99.81*table.n;i++)
1083 r = i*table.r/(100*table.n);
1085 expterm = exp(-0.25*r2);
1087 Vtab = 1/sqrt(r2+expterm);
1088 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1091 evaluate_table(table.tab,0,4,table.scale,r,&y0,&yp);
1092 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);
1094 abs_error_r=fabs(y0-Vtab);
1095 abs_error_r2=fabs(yp-(-1)*Ftab);
1097 rel_error_r=abs_error_r/y0;
1098 rel_error_r2=fabs(abs_error_r2/yp);
1101 if(rel_error_r>rel_error_r_old)
1103 rel_error_r_old=rel_error_r;
1107 if(rel_error_r2>rel_error_r2_old)
1109 rel_error_r2_old=rel_error_r2;
1114 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);
1115 printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1118 done_tabledata(&(td[0]));
1126 bondedtable_t make_bonded_table(FILE *fplog,char *fn,int angle)
1137 read_tables(fplog,fn,1,angle,&td);
1139 /* Convert the table from degrees to radians */
1140 for(i=0; i<td.nx; i++) {
1144 td.tabscale *= RAD2DEG;
1147 tab.scale = td.tabscale;
1148 snew(tab.tab,tab.n*4);
1149 copy2table(tab.n,0,4,td.x,td.v,td.f,tab.tab);
1150 done_tabledata(&td);