1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
45 #include "gmx_fatal.h"
56 /* All the possible (implemented) table functions */
81 /** Evaluates to true if the table type contains user data. */
82 #define ETAB_USER(e) ((e) == etabUSER || \
83 (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
90 /* This structure holds name and a flag that tells whether
91 this is a Coulomb type funtion */
92 static const t_tab_props tprops
[etabNR
] = {
95 { "LJ6Shift", FALSE
},
96 { "LJ12Shift", FALSE
},
102 { "Ewald-Switch", TRUE
},
103 { "Ewald-User", TRUE
},
104 { "Ewald-User-Switch", TRUE
},
105 { "LJ6Switch", FALSE
},
106 { "LJ12Switch", FALSE
},
107 { "COULSwitch", TRUE
},
108 { "LJ6-Encad shift", FALSE
},
109 { "LJ12-Encad shift", FALSE
},
110 { "COUL-Encad shift", TRUE
},
115 /* Index in the table that says which function to use */
116 enum { etiCOUL
, etiLJ6
, etiLJ12
, etiNR
};
124 #define pow2(x) ((x)*(x))
125 #define pow3(x) ((x)*(x)*(x))
126 #define pow4(x) ((x)*(x)*(x)*(x))
127 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
130 static double v_ewald_lr(double beta
,double r
)
134 return beta
*2/sqrt(M_PI
);
138 return gmx_erfd(beta
*r
)/r
;
142 void table_spline3_fill_ewald_lr(real
*tabf
,real
*tabv
,
143 int ntab
,int tableformat
,
150 gmx_bool bOutOfRange
;
151 double v_r0
,v_r1
,v_inrange
,vi
,a0
,a1
,a2dx
;
156 gmx_fatal(FARGS
,"Can not make a spline table with less than 2 points");
159 /* We need some margin to be able to divide table values by r
160 * in the kernel and also to do the integration arithmetics
161 * without going out of range. Furthemore, we divide by dx below.
163 tab_max
= GMX_REAL_MAX
*0.0001;
165 /* This function produces a table with:
166 * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
167 * maximum force error: V'''/(6*4)*dx^2
168 * The rms force error is the max error times 1/sqrt(5)=0.45.
173 case tableformatF
: stride
= 1; break;
174 case tableformatFDV0
: stride
= 4; break;
175 default: gmx_incons("Unknown table format");
182 for(i
=ntab
-1; i
>=0; i
--)
186 v_r0
= v_ewald_lr(beta
,x_r0
);
197 /* Linear continuation for the last point in range */
198 vi
= v_inrange
- dc
*(i
- i_inrange
)*dx
;
209 case tableformatFDV0
:
210 tabf
[i
*stride
+2] = vi
;
211 tabf
[i
*stride
+3] = 0;
214 gmx_incons("Unknown table format");
222 /* Get the potential at table point i-1 */
223 v_r1
= v_ewald_lr(beta
,(i
-1)*dx
);
225 if (v_r1
!= v_r1
|| v_r1
< -tab_max
|| v_r1
> tab_max
)
232 /* Calculate the average second derivative times dx over interval i-1 to i.
233 * Using the function values at the end points and in the middle.
235 a2dx
= (v_r0
+ v_r1
- 2*v_ewald_lr(beta
,x_r0
-0.5*dx
))/(0.25*dx
);
236 /* Set the derivative of the spline to match the difference in potential
237 * over the interval plus the average effect of the quadratic term.
238 * This is the essential step for minimizing the error in the force.
240 dc
= (v_r0
- v_r1
)/dx
+ 0.5*a2dx
;
245 /* Fill the table with the force, minus the derivative of the spline */
246 tabf
[i
*stride
] = -dc
;
250 /* tab[i] will contain the average of the splines over the two intervals */
251 tabf
[i
*stride
] += -0.5*dc
;
256 /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
257 * matching the potential at the two end points
258 * and the derivative dc at the end point xr.
262 a2dx
= (a1
*dx
+ v_r1
- a0
)*2/dx
;
264 /* Set dc to the derivative at the next point */
267 if (dc_new
!= dc_new
|| dc_new
< -tab_max
|| dc_new
> tab_max
)
277 tabf
[(i
-1)*stride
] = -0.5*dc
;
279 /* Currently the last value only contains half the force: double it */
282 if (tableformat
== tableformatFDV0
)
284 /* Store the force difference in the second entry */
285 for(i
=0; i
<ntab
-1; i
++)
287 tabf
[i
*stride
+1] = tabf
[(i
+1)*stride
] - tabf
[i
*stride
];
289 tabf
[(ntab
-1)*stride
+1] = -tabf
[i
*stride
];
293 /* The scale (1/spacing) for third order spline interpolation
294 * of the Ewald mesh contribution which needs to be subtracted
295 * from the non-bonded interactions.
297 real
ewald_spline3_table_scale(real ewaldcoeff
,real rc
)
299 double erf_x_d3
=1.0522; /* max of (erf(x)/x)''' */
303 /* Force tolerance: single precision accuracy */
304 ftol
= GMX_FLOAT_EPS
;
305 sc_f
= sqrt(erf_x_d3
/(6*4*ftol
*ewaldcoeff
))*ewaldcoeff
;
307 /* Energy tolerance: 10x more accurate than the cut-off jump */
308 etol
= 0.1*gmx_erfc(ewaldcoeff
*rc
);
309 etol
= max(etol
,GMX_REAL_EPS
);
310 sc_e
= pow(erf_x_d3
/(6*12*sqrt(3)*etol
),1.0/3.0)*ewaldcoeff
;
312 return max(sc_f
,sc_e
);
315 /* Calculate the potential and force for an r value
316 * in exactly the same way it is done in the inner loop.
317 * VFtab is a pointer to the table data, offset is
318 * the point where we should begin and stride is
319 * 4 if we have a buckingham table, 3 otherwise.
320 * If you want to evaluate table no N, set offset to 4*N.
322 * We use normal precision here, since that is what we
323 * will use in the inner loops.
325 static void evaluate_table(real VFtab
[], int offset
, int stride
,
326 real tabscale
, real r
, real
*y
, real
*yp
)
330 real Y
,F
,Geps
,Heps2
,Fp
;
339 Geps
= eps
*VFtab
[n
+2];
340 Heps2
= eps2
*VFtab
[n
+3];
343 *yp
= (Fp
+Geps
+2.0*Heps2
)*tabscale
;
346 static void copy2table(int n
,int offset
,int stride
,
347 double x
[],double Vtab
[],double Ftab
[],
350 /* Use double prec. for the intermediary variables
351 * and temporary x/vtab/vtab2 data to avoid unnecessary
358 for(i
=0; (i
<n
); i
++) {
362 G
= 3*(Vtab
[i
+1] - Vtab
[i
]) + (Ftab
[i
+1] + 2*Ftab
[i
])*h
;
363 H
= -2*(Vtab
[i
+1] - Vtab
[i
]) - (Ftab
[i
+1] + Ftab
[i
])*h
;
365 /* Fill the last entry with a linear potential,
366 * this is mainly for rounding issues with angle and dihedral potentials.
372 nn0
= offset
+ i
*stride
;
380 static void init_table(FILE *fp
,int n
,int nx0
,
381 double tabscale
,t_tabledata
*td
,gmx_bool bAlloc
)
387 td
->tabscale
= tabscale
;
393 for(i
=0; (i
<td
->nx
); i
++)
394 td
->x
[i
] = i
/tabscale
;
397 static void spline_forces(int nx
,double h
,double v
[],gmx_bool bS3
,gmx_bool bE3
,
404 /* Formulas can be found in:
405 * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
408 if (nx
< 4 && (bS3
|| bE3
))
409 gmx_fatal(FARGS
,"Can not generate splines with third derivative boundary conditions with less than 4 (%d) points",nx
);
411 /* To make life easy we initially set the spacing to 1
412 * and correct for this at the end.
416 /* Fit V''' at the start */
417 v3
= v
[3] - 3*v
[2] + 3*v
[1] - v
[0];
419 fprintf(debug
,"The left third derivative is %g\n",v3
/(h
*h
*h
));
420 b_s
= 2*(v
[1] - v
[0]) + v3
/6;
424 /* Fit V'' at the start */
427 v2
= -v
[3] + 4*v
[2] - 5*v
[1] + 2*v
[0];
428 /* v2 = v[2] - 2*v[1] + v[0]; */
430 fprintf(debug
,"The left second derivative is %g\n",v2
/(h
*h
));
431 b_s
= 3*(v
[1] - v
[0]) - v2
/2;
435 b_s
= 3*(v
[2] - v
[0]) + f
[0]*h
;
439 /* Fit V''' at the end */
440 v3
= v
[nx
-1] - 3*v
[nx
-2] + 3*v
[nx
-3] - v
[nx
-4];
442 fprintf(debug
,"The right third derivative is %g\n",v3
/(h
*h
*h
));
443 b_e
= 2*(v
[nx
-1] - v
[nx
-2]) + v3
/6;
446 /* V'=0 at the end */
447 b_e
= 3*(v
[nx
-1] - v
[nx
-3]) + f
[nx
-1]*h
;
452 beta
= (bS3
? 1 : 4);
454 /* For V'' fitting */
455 /* beta = (bS3 ? 2 : 4); */
458 for(i
=start
+1; i
<end
; i
++) {
461 b
= 3*(v
[i
+1] - v
[i
-1]);
462 f
[i
] = (b
- f
[i
-1])/beta
;
464 gamma
[end
-1] = 1/beta
;
465 beta
= (bE3
? 1 : 4) - gamma
[end
-1];
466 f
[end
-1] = (b_e
- f
[end
-2])/beta
;
468 for(i
=end
-2; i
>=start
; i
--)
469 f
[i
] -= gamma
[i
+1]*f
[i
+1];
472 /* Correct for the minus sign and the spacing */
473 for(i
=start
; i
<end
; i
++)
477 static void set_forces(FILE *fp
,int angle
,
478 int nx
,double h
,double v
[],double f
[],
485 "Force generation for dihedral tables is not (yet) implemented");
488 while (v
[start
] == 0)
500 fprintf(fp
,"Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
501 table
+1,start
*h
,end
==nx
? "V'''" : "V'=0",(end
-1)*h
);
502 spline_forces(end
-start
,h
,v
+start
,TRUE
,end
==nx
,f
+start
);
505 static void read_tables(FILE *fp
,const char *fn
,
506 int ntab
,int angle
,t_tabledata td
[])
510 double **yy
=NULL
,start
,end
,dx0
,dx1
,ssd
,vm
,vp
,f
,numf
;
511 int k
,i
,nx
,nx0
=0,ny
,nny
,ns
;
512 gmx_bool bAllZero
,bZeroV
,bZeroF
;
516 libfn
= gmxlibfn(fn
);
517 nx
= read_xvg(libfn
,&yy
,&ny
);
519 gmx_fatal(FARGS
,"Trying to read file %s, but nr columns = %d, should be %d",
524 "The first distance in file %s is %f nm instead of %f nm",
532 if (yy
[0][0] != start
|| yy
[0][nx
-1] != end
)
533 gmx_fatal(FARGS
,"The angles in file %s should go from %f to %f instead of %f to %f\n",
534 libfn
,start
,end
,yy
[0][0],yy
[0][nx
-1]);
537 tabscale
= (nx
-1)/(yy
[0][nx
-1] - yy
[0][0]);
540 fprintf(fp
,"Read user tables from %s with %d data points.\n",libfn
,nx
);
542 fprintf(fp
,"Tabscale = %g points/nm\n",tabscale
);
546 for(k
=0; k
<ntab
; k
++) {
549 for(i
=0; (i
< nx
); i
++) {
551 dx0
= yy
[0][i
-1] - yy
[0][i
-2];
552 dx1
= yy
[0][i
] - yy
[0][i
-1];
553 /* Check for 1% deviation in spacing */
554 if (fabs(dx1
- dx0
) >= 0.005*(fabs(dx0
) + fabs(dx1
))) {
555 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
]);
558 if (yy
[1+k
*2][i
] != 0) {
564 if (yy
[1+k
*2][i
] > 0.01*GMX_REAL_MAX
||
565 yy
[1+k
*2][i
] < -0.01*GMX_REAL_MAX
) {
566 gmx_fatal(FARGS
,"Out of range potential value %g in file '%s'",
570 if (yy
[1+k
*2+1][i
] != 0) {
576 if (yy
[1+k
*2+1][i
] > 0.01*GMX_REAL_MAX
||
577 yy
[1+k
*2+1][i
] < -0.01*GMX_REAL_MAX
) {
578 gmx_fatal(FARGS
,"Out of range force value %g in file '%s'",
584 if (!bZeroV
&& bZeroF
) {
585 set_forces(fp
,angle
,nx
,1/tabscale
,yy
[1+k
*2],yy
[1+k
*2+1],k
);
587 /* Check if the second column is close to minus the numerical
588 * derivative of the first column.
592 for(i
=1; (i
< nx
-1); i
++) {
596 if (vm
!= 0 && vp
!= 0 && f
!= 0) {
597 /* Take the centered difference */
598 numf
= -(vp
- vm
)*0.5*tabscale
;
599 ssd
+= fabs(2*(f
- numf
)/(f
+ numf
));
605 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));
607 fprintf(debug
,"%s",buf
);
610 fprintf(fp
,"\nWARNING: %s\n",buf
);
611 fprintf(stderr
,"\nWARNING: %s\n",buf
);
616 if (bAllZero
&& fp
) {
617 fprintf(fp
,"\nNOTE: All elements in table %s are zero\n\n",libfn
);
620 for(k
=0; (k
<ntab
); k
++) {
621 init_table(fp
,nx
,nx0
,tabscale
,&(td
[k
]),TRUE
);
622 for(i
=0; (i
<nx
); i
++) {
623 td
[k
].x
[i
] = yy
[0][i
];
624 td
[k
].v
[i
] = yy
[2*k
+1][i
];
625 td
[k
].f
[i
] = yy
[2*k
+2][i
];
628 for(i
=0; (i
<ny
); i
++)
634 static void done_tabledata(t_tabledata
*td
)
646 static void fill_table(t_tabledata
*td
,int tp
,const t_forcerec
*fr
)
648 /* Fill the table according to the formulas in the manual.
649 * In principle, we only need the potential and the second
650 * derivative, but then we would have to do lots of calculations
651 * in the inner loop. By precalculating some terms (see manual)
652 * we get better eventual performance, despite a larger table.
654 * Since some of these higher-order terms are very small,
655 * we always use double precision to calculate them here, in order
656 * to avoid unnecessary loss of precision.
663 double r1
,rc
,r12
,r13
;
665 double expr
,Vtab
,Ftab
;
666 /* Parameters for David's function */
667 double A
=0,B
=0,C
=0,A_3
=0,B_4
=0;
668 /* Parameters for the switching function */
670 /* Temporary parameters */
671 gmx_bool bSwitch
,bShift
;
672 double ewc
=fr
->ewaldcoeff
;
673 double isp
= 0.564189583547756;
675 bSwitch
= ((tp
== etabLJ6Switch
) || (tp
== etabLJ12Switch
) ||
676 (tp
== etabCOULSwitch
) ||
677 (tp
== etabEwaldSwitch
) || (tp
== etabEwaldUserSwitch
));
678 bShift
= ((tp
== etabLJ6Shift
) || (tp
== etabLJ12Shift
) ||
683 if (tprops
[tp
].bCoulomb
) {
684 r1
= fr
->rcoulomb_switch
;
688 r1
= fr
->rvdw_switch
;
692 ksw
= 1.0/(pow5(rc
-r1
));
698 else if (tp
== etabLJ6Shift
)
703 A
= p
* ((p
+1)*r1
-(p
+4)*rc
)/(pow(rc
,p
+2)*pow2(rc
-r1
));
704 B
= -p
* ((p
+1)*r1
-(p
+3)*rc
)/(pow(rc
,p
+2)*pow3(rc
-r1
));
705 C
= 1.0/pow(rc
,p
)-A
/3.0*pow3(rc
-r1
)-B
/4.0*pow4(rc
-r1
);
706 if (tp
== etabLJ6Shift
) {
714 if (debug
) { fprintf(debug
,"Setting up tables\n"); fflush(debug
); }
717 fp
=xvgropen("switch.xvg","switch","r","s");
720 for(i
=td
->nx0
; (i
<td
->nx
); i
++) {
724 if (gmx_within_tol(reppow
,12.0,10*GMX_DOUBLE_EPS
)) {
727 r12
= pow(r
,-reppow
);
732 /* swi is function, swi1 1st derivative and swi2 2nd derivative */
733 /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
734 * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
736 * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
745 swi
= 1 - 10*pow3(r
-r1
)*ksw
*pow2(rc
-r1
)
746 + 15*pow4(r
-r1
)*ksw
*(rc
-r1
) - 6*pow5(r
-r1
)*ksw
;
747 swi1
= -30*pow2(r
-r1
)*ksw
*pow2(rc
-r1
)
748 + 60*pow3(r
-r1
)*ksw
*(rc
-r1
) - 30*pow4(r
-r1
)*ksw
;
751 else { /* not really needed, but avoids compiler warnings... */
756 fprintf(fp
,"%10g %10g %10g %10g\n",r
,swi
,swi1
,swi2
);
779 Ftab
= reppow
*Vtab
/r
;
786 Ftab
= reppow
*Vtab
/r
;
791 Vtab
= -(r6
-6.0*(rc
-r
)*rc6
/rc
-rc6
);
792 Ftab
= -(6.0*r6
/r
-6.0*rc6
/rc
);
800 Vtab
= r12
-12.0*(rc
-r
)*rc6
*rc6
/rc
-1.0*rc6
*rc6
;
801 Ftab
= 12.0*r12
/r
-12.0*rc6
*rc6
/rc
;
819 case etabEwaldSwitch
:
820 Vtab
= gmx_erfc(ewc
*r
)/r
;
821 Ftab
= gmx_erfc(ewc
*r
)/r2
+2*exp(-(ewc
*ewc
*r2
))*ewc
*isp
/r
;
824 case etabEwaldUserSwitch
:
825 /* Only calculate minus the reciprocal space contribution */
826 Vtab
= -gmx_erf(ewc
*r
)/r
;
827 Ftab
= -gmx_erf(ewc
*r
)/r2
+2*exp(-(ewc
*ewc
*r2
))*ewc
*isp
/r
;
831 Vtab
= 1.0/r
+ fr
->k_rf
*r2
- fr
->c_rf
;
832 Ftab
= 1.0/r2
- 2*fr
->k_rf
*r
;
833 if (tp
== etabRF_ZERO
&& r
>= rc
) {
845 Vtab
= 1.0/r
-(rc
-r
)/(rc
*rc
)-1.0/rc
;
846 Ftab
= 1.0/r2
-1.0/(rc
*rc
);
853 gmx_fatal(FARGS
,"Table type %d not implemented yet. (%s,%d)",
854 tp
,__FILE__
,__LINE__
);
857 /* Normal coulomb with cut-off correction for potential */
860 /* If in Shifting range add something to it */
864 Vtab
+= - A_3
*r13
- B_4
*r12
*r12
;
865 Ftab
+= A
*r12
+ B
*r13
;
875 if ((r
> r1
) && bSwitch
) {
876 Ftab
= Ftab
*swi
- Vtab
*swi1
;
880 /* Convert to single precision when we store to mem */
885 /* Continue the table linearly from nx0 to 0.
886 * These values are only required for energy minimization with overlap or TPI.
888 for(i
=td
->nx0
-1; i
>=0; i
--) {
889 td
->v
[i
] = td
->v
[i
+1] + td
->f
[i
+1]*(td
->x
[i
+1] - td
->x
[i
]);
890 td
->f
[i
] = td
->f
[i
+1];
898 static void set_table_type(int tabsel
[],const t_forcerec
*fr
,gmx_bool b14only
)
902 /* Set the different table indices.
908 switch (fr
->eeltype
) {
914 case eelPMEUSERSWITCH
:
921 eltype
= fr
->eeltype
;
926 tabsel
[etiCOUL
] = etabCOUL
;
929 tabsel
[etiCOUL
] = etabShift
;
932 if (fr
->rcoulomb
> fr
->rcoulomb_switch
)
933 tabsel
[etiCOUL
] = etabShift
;
935 tabsel
[etiCOUL
] = etabCOUL
;
940 tabsel
[etiCOUL
] = etabEwald
;
943 tabsel
[etiCOUL
] = etabEwaldSwitch
;
946 tabsel
[etiCOUL
] = etabEwaldUser
;
948 case eelPMEUSERSWITCH
:
949 tabsel
[etiCOUL
] = etabEwaldUserSwitch
;
954 tabsel
[etiCOUL
] = etabRF
;
957 tabsel
[etiCOUL
] = etabRF_ZERO
;
960 tabsel
[etiCOUL
] = etabCOULSwitch
;
963 tabsel
[etiCOUL
] = etabUSER
;
966 tabsel
[etiCOUL
] = etabCOULEncad
;
969 gmx_fatal(FARGS
,"Invalid eeltype %d",eltype
);
972 /* Van der Waals time */
973 if (fr
->bBHAM
&& !b14only
) {
974 tabsel
[etiLJ6
] = etabLJ6
;
975 tabsel
[etiLJ12
] = etabEXPMIN
;
977 if (b14only
&& fr
->vdwtype
!= evdwUSER
)
980 vdwtype
= fr
->vdwtype
;
984 tabsel
[etiLJ6
] = etabLJ6Switch
;
985 tabsel
[etiLJ12
] = etabLJ12Switch
;
988 tabsel
[etiLJ6
] = etabLJ6Shift
;
989 tabsel
[etiLJ12
] = etabLJ12Shift
;
992 tabsel
[etiLJ6
] = etabUSER
;
993 tabsel
[etiLJ12
] = etabUSER
;
996 tabsel
[etiLJ6
] = etabLJ6
;
997 tabsel
[etiLJ12
] = etabLJ12
;
1000 tabsel
[etiLJ6
] = etabLJ6Encad
;
1001 tabsel
[etiLJ12
] = etabLJ12Encad
;
1004 gmx_fatal(FARGS
,"Invalid vdwtype %d in %s line %d",vdwtype
,
1010 t_forcetable
make_tables(FILE *out
,const output_env_t oenv
,
1011 const t_forcerec
*fr
,
1012 gmx_bool bVerbose
,const char *fn
,
1013 real rtab
,int flags
)
1015 const char *fns
[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
1016 const char *fns14
[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" };
1019 gmx_bool b14only
,bReadTab
,bGenTab
;
1021 int i
,j
,k
,nx
,nx0
,tabsel
[etiNR
];
1025 b14only
= (flags
& GMX_MAKETABLES_14ONLY
);
1027 if (flags
& GMX_MAKETABLES_FORCEUSER
) {
1028 tabsel
[etiCOUL
] = etabUSER
;
1029 tabsel
[etiLJ6
] = etabUSER
;
1030 tabsel
[etiLJ12
] = etabUSER
;
1032 set_table_type(tabsel
,fr
,b14only
);
1038 table
.scale_exp
= 0;
1042 /* Check whether we have to read or generate */
1045 for(i
=0; (i
<etiNR
); i
++) {
1046 if (ETAB_USER(tabsel
[i
]))
1048 if (tabsel
[i
] != etabUSER
)
1052 read_tables(out
,fn
,etiNR
,0,td
);
1053 if (rtab
== 0 || (flags
& GMX_MAKETABLES_14ONLY
)) {
1054 rtab
= td
[0].x
[td
[0].nx
-1];
1058 if (td
[0].x
[td
[0].nx
-1] < rtab
)
1059 gmx_fatal(FARGS
,"Tables in file %s not long enough for cut-off:\n"
1060 "\tshould be at least %f nm\n",fn
,rtab
);
1061 nx
= table
.n
= (int)(rtab
*td
[0].tabscale
+ 0.5);
1063 table
.scale
= td
[0].tabscale
;
1069 table
.scale
= 2000.0;
1071 table
.scale
= 500.0;
1073 nx
= table
.n
= rtab
*table
.scale
;
1077 if(fr
->bham_b_max
!=0)
1078 table
.scale_exp
= table
.scale
/fr
->bham_b_max
;
1080 table
.scale_exp
= table
.scale
;
1083 /* Each table type (e.g. coul,lj6,lj12) requires four
1084 * numbers per nx+1 data points. For performance reasons we want
1085 * the table data to be aligned to 16-byte.
1087 snew_aligned(table
.tab
, 12*(nx
+1)*sizeof(real
),16);
1089 for(k
=0; (k
<etiNR
); k
++) {
1090 if (tabsel
[k
] != etabUSER
) {
1091 init_table(out
,nx
,nx0
,
1092 (tabsel
[k
] == etabEXPMIN
) ? table
.scale_exp
: table
.scale
,
1093 &(td
[k
]),!bReadTab
);
1094 fill_table(&(td
[k
]),tabsel
[k
],fr
);
1096 fprintf(out
,"%s table with %d data points for %s%s.\n"
1097 "Tabscale = %g points/nm\n",
1098 ETAB_USER(tabsel
[k
]) ? "Modified" : "Generated",
1099 td
[k
].nx
,b14only
?"1-4 ":"",tprops
[tabsel
[k
]].name
,
1102 copy2table(table
.n
,k
*4,12,td
[k
].x
,td
[k
].v
,td
[k
].f
,table
.tab
);
1104 if (bDebugMode() && bVerbose
) {
1106 fp
=xvgropen(fns14
[k
],fns14
[k
],"r","V",oenv
);
1108 fp
=xvgropen(fns
[k
],fns
[k
],"r","V",oenv
);
1109 /* plot the output 5 times denser than the table data */
1110 for(i
=5*((nx0
+1)/2); i
<5*table
.n
; i
++) {
1111 x0
= i
*table
.r
/(5*(table
.n
-1));
1112 evaluate_table(table
.tab
,4*k
,12,table
.scale
,x0
,&y0
,&yp
);
1113 fprintf(fp
,"%15.10e %15.10e %15.10e\n",x0
,y0
,yp
);
1117 done_tabledata(&(td
[k
]));
1124 t_forcetable
make_gb_table(FILE *out
,const output_env_t oenv
,
1125 const t_forcerec
*fr
,
1129 const char *fns
[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
1130 const char *fns14
[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
1133 gmx_bool bReadTab
,bGenTab
;
1135 int i
,j
,k
,nx
,nx0
,tabsel
[etiNR
];
1136 double r
,r2
,Vtab
,Ftab
,expterm
;
1140 double abs_error_r
, abs_error_r2
;
1141 double rel_error_r
, rel_error_r2
;
1142 double rel_error_r_old
=0, rel_error_r2_old
=0;
1143 double x0_r_error
, x0_r2_error
;
1146 /* Only set a Coulomb table for GB */
1153 /* Set the table dimensions for GB, not really necessary to
1154 * use etiNR (since we only have one table, but ...)
1157 table
.r
= fr
->gbtabr
;
1158 table
.scale
= fr
->gbtabscale
;
1159 table
.scale_exp
= 0;
1160 table
.n
= table
.scale
*table
.r
;
1162 nx
= table
.scale
*table
.r
;
1164 /* Check whether we have to read or generate
1165 * We will always generate a table, so remove the read code
1166 * (Compare with original make_table function
1171 /* Each table type (e.g. coul,lj6,lj12) requires four
1172 * numbers per datapoint. For performance reasons we want
1173 * the table data to be aligned to 16-byte. This is accomplished
1174 * by allocating 16 bytes extra to a temporary pointer, and then
1175 * calculating an aligned pointer. This new pointer must not be
1176 * used in a free() call, but thankfully we're sloppy enough not
1180 snew_aligned(table
.tab
,4*nx
,16);
1182 init_table(out
,nx
,nx0
,table
.scale
,&(td
[0]),!bReadTab
);
1184 /* Local implementation so we don't have to use the etabGB
1185 * enum above, which will cause problems later when
1186 * making the other tables (right now even though we are using
1187 * GB, the normal Coulomb tables will be created, but this
1188 * will cause a problem since fr->eeltype==etabGB which will not
1189 * be defined in fill_table and set_table_type
1198 expterm
= exp(-0.25*r2
);
1200 Vtab
= 1/sqrt(r2
+expterm
);
1201 Ftab
= (r
-0.25*r
*expterm
)/((r2
+expterm
)*sqrt(r2
+expterm
));
1203 /* Convert to single precision when we store to mem */
1209 copy2table(table
.n
,0,4,td
[0].x
,td
[0].v
,td
[0].f
,table
.tab
);
1213 fp
=xvgropen(fns
[0],fns
[0],"r","V",oenv
);
1214 /* plot the output 5 times denser than the table data */
1215 /* for(i=5*nx0;i<5*table.n;i++) */
1216 for(i
=nx0
;i
<table
.n
;i
++)
1218 /* x0=i*table.r/(5*table.n); */
1219 x0
=i
*table
.r
/table
.n
;
1220 evaluate_table(table
.tab
,0,4,table
.scale
,x0
,&y0
,&yp
);
1221 fprintf(fp
,"%15.10e %15.10e %15.10e\n",x0
,y0
,yp
);
1228 for(i=100*nx0;i<99.81*table.n;i++)
1230 r = i*table.r/(100*table.n);
1232 expterm = exp(-0.25*r2);
1234 Vtab = 1/sqrt(r2+expterm);
1235 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1238 evaluate_table(table.tab,0,4,table.scale,r,&y0,&yp);
1239 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);
1241 abs_error_r=fabs(y0-Vtab);
1242 abs_error_r2=fabs(yp-(-1)*Ftab);
1244 rel_error_r=abs_error_r/y0;
1245 rel_error_r2=fabs(abs_error_r2/yp);
1248 if(rel_error_r>rel_error_r_old)
1250 rel_error_r_old=rel_error_r;
1254 if(rel_error_r2>rel_error_r2_old)
1256 rel_error_r2_old=rel_error_r2;
1261 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);
1262 printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1265 done_tabledata(&(td
[0]));
1273 t_forcetable
make_atf_table(FILE *out
,const output_env_t oenv
,
1274 const t_forcerec
*fr
,
1278 const char *fns
[3] = { "tf_tab.xvg", "atfdtab.xvg", "atfrtab.xvg" };
1283 real rx
, ry
, rz
, box_r
;
1288 /* Set the table dimensions for ATF, not really necessary to
1289 * use etiNR (since we only have one table, but ...)
1293 if (fr
->adress_type
== eAdressSphere
){
1294 /* take half box diagonal direction as tab range */
1295 rx
= 0.5*box
[0][0]+0.5*box
[1][0]+0.5*box
[2][0];
1296 ry
= 0.5*box
[0][1]+0.5*box
[1][1]+0.5*box
[2][1];
1297 rz
= 0.5*box
[0][2]+0.5*box
[1][2]+0.5*box
[2][2];
1298 box_r
= sqrt(rx
*rx
+ry
*ry
+rz
*rz
);
1301 /* xsplit: take half box x direction as tab range */
1302 box_r
= box
[0][0]/2;
1307 table
.scale_exp
= 0;
1311 read_tables(out
,fn
,1,0,td
);
1312 rtab
= td
[0].x
[td
[0].nx
-1];
1314 if (fr
->adress_type
== eAdressXSplit
&& (rtab
< box
[0][0]/2)){
1315 gmx_fatal(FARGS
,"AdResS full box therm force table in file %s extends to %f:\n"
1316 "\tshould extend to at least half the length of the box in x-direction"
1317 "%f\n",fn
,rtab
, box
[0][0]/2);
1320 gmx_fatal(FARGS
,"AdResS full box therm force table in file %s extends to %f:\n"
1321 "\tshould extend to at least for spherical adress"
1322 "%f (=distance from center to furthermost point in box \n",fn
,rtab
, box_r
);
1328 table
.scale
= td
[0].tabscale
;
1331 /* Each table type (e.g. coul,lj6,lj12) requires four
1332 * numbers per datapoint. For performance reasons we want
1333 * the table data to be aligned to 16-byte. This is accomplished
1334 * by allocating 16 bytes extra to a temporary pointer, and then
1335 * calculating an aligned pointer. This new pointer must not be
1336 * used in a free() call, but thankfully we're sloppy enough not
1340 snew_aligned(table
.tab
,4*nx
,16);
1342 copy2table(table
.n
,0,4,td
[0].x
,td
[0].v
,td
[0].f
,table
.tab
);
1346 fp
=xvgropen(fns
[0],fns
[0],"r","V",oenv
);
1347 /* plot the output 5 times denser than the table data */
1348 /* for(i=5*nx0;i<5*table.n;i++) */
1350 for(i
=5*((nx0
+1)/2); i
<5*table
.n
; i
++)
1352 /* x0=i*table.r/(5*table.n); */
1353 x0
= i
*table
.r
/(5*(table
.n
-1));
1354 evaluate_table(table
.tab
,0,4,table
.scale
,x0
,&y0
,&yp
);
1355 fprintf(fp
,"%15.10e %15.10e %15.10e\n",x0
,y0
,yp
);
1361 done_tabledata(&(td
[0]));
1367 bondedtable_t
make_bonded_table(FILE *fplog
,char *fn
,int angle
)
1378 read_tables(fplog
,fn
,1,angle
,&td
);
1380 /* Convert the table from degrees to radians */
1381 for(i
=0; i
<td
.nx
; i
++) {
1385 td
.tabscale
*= RAD2DEG
;
1388 tab
.scale
= td
.tabscale
;
1389 snew(tab
.tab
,tab
.n
*4);
1390 copy2table(tab
.n
,0,4,td
.x
,td
.v
,td
.f
,tab
.tab
);
1391 done_tabledata(&td
);