}
}
-void table_spline3_fill_ewald_lr(real *tabf,real *tabv,
- int ntab,int tableformat,
- real dx,real beta)
+void table_spline3_fill_ewald_lr(real *table_f,
+ real *table_v,
+ real *table_fdv0,
+ int ntab,
+ real dx,
+ real beta)
{
real tab_max;
- int stride=0;
int i,i_inrange;
double dc,dc_new;
gmx_bool bOutOfRange;
* The rms force error is the max error times 1/sqrt(5)=0.45.
*/
- switch (tableformat)
- {
- case tableformatF: stride = 1; break;
- case tableformatFDV0: stride = 4; break;
- default: gmx_incons("Unknown table format");
- }
-
bOutOfRange = FALSE;
i_inrange = ntab;
v_inrange = 0;
vi = v_inrange - dc*(i - i_inrange)*dx;
}
- switch (tableformat)
+ if(table_v!=NULL)
{
- case tableformatF:
- if (tabv != NULL)
- {
- tabv[i] = vi;
- }
- break;
- case tableformatFDV0:
- tabf[i*stride+2] = vi;
- tabf[i*stride+3] = 0;
- break;
- default:
- gmx_incons("Unknown table format");
+ table_v[i] = vi;
}
if (i == 0)
if (i == ntab - 1)
{
/* Fill the table with the force, minus the derivative of the spline */
- tabf[i*stride] = -dc;
+ table_f[i] = -dc;
}
else
{
/* tab[i] will contain the average of the splines over the two intervals */
- tabf[i*stride] += -0.5*dc;
+ table_f[i] += -0.5*dc;
}
if (!bOutOfRange)
}
}
- tabf[(i-1)*stride] = -0.5*dc;
+ table_f[(i-1)] = -0.5*dc;
}
/* Currently the last value only contains half the force: double it */
- tabf[0] *= 2;
+ table_f[0] *= 2;
- if (tableformat == tableformatFDV0)
+ if(table_v!=NULL && table_fdv0!=NULL)
{
- /* Store the force difference in the second entry */
- for(i=0; i<ntab-1; i++)
+ /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
+ * init_ewald_f_table().
+ */
+ for(i=0;i<ntab-1;i++)
{
- tabf[i*stride+1] = tabf[(i+1)*stride] - tabf[i*stride];
+ table_fdv0[4*i] = table_f[i];
+ table_fdv0[4*i+1] = table_f[i+1]-table_f[i];
+ table_fdv0[4*i+2] = table_v[i];
+ table_fdv0[4*i+3] = 0.0;
}
- tabf[(ntab-1)*stride+1] = -tabf[i*stride];
+ table_fdv0[4*(ntab-1)] = table_f[(ntab-1)];
+ table_fdv0[4*(ntab-1)+1] = -table_f[(ntab-1)];
+ table_fdv0[4*(ntab-1)+2] = table_v[(ntab-1)];
+ table_fdv0[4*(ntab-1)+3] = 0.0;
}
}
}
static void copy2table(int n,int offset,int stride,
- double x[],double Vtab[],double Ftab[],
+ double x[],double Vtab[],double Ftab[],real scalefactor,
real dest[])
{
/* Use double prec. for the intermediary variables
H = 0;
}
nn0 = offset + i*stride;
- dest[nn0] = Vtab[i];
- dest[nn0+1] = F;
- dest[nn0+2] = G;
- dest[nn0+3] = H;
+ dest[nn0] = scalefactor*Vtab[i];
+ dest[nn0+1] = scalefactor*F;
+ dest[nn0+2] = scalefactor*G;
+ dest[nn0+3] = scalefactor*H;
}
}
switch (tp) {
case etabLJ6:
- /* Dispersion */
- Vtab = -r6;
- Ftab = 6.0*Vtab/r;
- break;
+ /* Dispersion */
+ Vtab = -r6;
+ Ftab = 6.0*Vtab/r;
+ break;
case etabLJ6Switch:
case etabLJ6Shift:
/* Dispersion */
if (r < rc) {
- Vtab = -r6;
- Ftab = 6.0*Vtab/r;
+ Vtab = -r6;
+ Ftab = 6.0*Vtab/r;
+ break;
}
break;
case etabLJ12:
- /* Repulsion */
- Vtab = r12;
- Ftab = reppow*Vtab/r;
+ /* Repulsion */
+ Vtab = r12;
+ Ftab = reppow*Vtab/r;
break;
case etabLJ12Switch:
case etabLJ12Shift:
/* Repulsion */
if (r < rc) {
- Vtab = r12;
- Ftab = reppow*Vtab/r;
- }
+ Vtab = r12;
+ Ftab = reppow*Vtab/r;
+ }
break;
case etabLJ6Encad:
if(r < rc) {
break;
case etabLJ12Encad:
if(r < rc) {
- Vtab = r12-12.0*(rc-r)*rc6*rc6/rc-1.0*rc6*rc6;
- Ftab = 12.0*r12/r-12.0*rc6*rc6/rc;
- } else { /* r>rc */
+ Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
+ Ftab = -(6.0*r6/r-6.0*rc6/rc);
+ } else { /* r>rc */
Vtab = 0;
Ftab = 0;
}
if ((r > r1) && bSwitch) {
Ftab = Ftab*swi - Vtab*swi1;
Vtab = Vtab*swi;
- }
-
+ }
+
/* Convert to single precision when we store to mem */
td->v[i] = Vtab;
td->f[i] = Ftab;
gmx_bool b14only,bReadTab,bGenTab;
real x0,y0,yp;
int i,j,k,nx,nx0,tabsel[etiNR];
-
+ real scalefactor;
+
t_forcetable table;
b14only = (flags & GMX_MAKETABLES_14ONLY);
nx0 = 10;
nx = 0;
+ table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
+ table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
+ table.formatsize = 4;
+ table.ninteractions = 3;
+ table.stride = table.formatsize*table.ninteractions;
+
/* Check whether we have to read or generate */
bReadTab = FALSE;
bGenTab = FALSE;
* numbers per nx+1 data points. For performance reasons we want
* the table data to be aligned to 16-byte.
*/
- snew_aligned(table.tab, 12*(nx+1)*sizeof(real),16);
+ snew_aligned(table.data, 12*(nx+1)*sizeof(real),16);
for(k=0; (k<etiNR); k++) {
if (tabsel[k] != etabUSER) {
td[k].nx,b14only?"1-4 ":"",tprops[tabsel[k]].name,
td[k].tabscale);
}
- copy2table(table.n,k*4,12,td[k].x,td[k].v,td[k].f,table.tab);
+
+ /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
+ * by including the derivative constants (6.0 or 12.0) in the parameters, since
+ * we no longer calculate force in most steps. This means the c6/c12 parameters
+ * have been scaled up, so we need to scale down the table interactions too.
+ * It comes here since we need to scale user tables too.
+ */
+ if(k==etiLJ6)
+ {
+ scalefactor = 1.0/6.0;
+ }
+ else if(k==etiLJ12 && tabsel[k]!=etabEXPMIN)
+ {
+ scalefactor = 1.0/12.0;
+ }
+ else
+ {
+ scalefactor = 1.0;
+ }
+
+ copy2table(table.n,k*4,12,td[k].x,td[k].v,td[k].f,scalefactor,table.data);
if (bDebugMode() && bVerbose) {
if (b14only)
/* plot the output 5 times denser than the table data */
for(i=5*((nx0+1)/2); i<5*table.n; i++) {
x0 = i*table.r/(5*(table.n-1));
- evaluate_table(table.tab,4*k,12,table.scale,x0,&y0,&yp);
+ evaluate_table(table.data,4*k,12,table.scale,x0,&y0,&yp);
fprintf(fp,"%15.10e %15.10e %15.10e\n",x0,y0,yp);
}
gmx_fio_fclose(fp);
* use etiNR (since we only have one table, but ...)
*/
snew(td,1);
- table.r = fr->gbtabr;
- table.scale = fr->gbtabscale;
- table.scale_exp = 0;
- table.n = table.scale*table.r;
- nx0 = 0;
- nx = table.scale*table.r;
+ table.interaction = GMX_TABLE_INTERACTION_ELEC;
+ table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
+ table.r = fr->gbtabr;
+ table.scale = fr->gbtabscale;
+ table.scale_exp = 0;
+ table.n = table.scale*table.r;
+ table.formatsize = 4;
+ table.ninteractions = 1;
+ table.stride = table.formatsize*table.ninteractions;
+ nx0 = 0;
+ nx = table.scale*table.r;
/* Check whether we have to read or generate
* We will always generate a table, so remove the read code
* to do this :-)
*/
- snew_aligned(table.tab,4*nx,16);
+ snew_aligned(table.data,4*nx,16);
init_table(out,nx,nx0,table.scale,&(td[0]),!bReadTab);
}
- copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,table.tab);
+ copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,1.0,table.data);
if(bDebugMode())
{
{
/* x0=i*table.r/(5*table.n); */
x0=i*table.r/table.n;
- evaluate_table(table.tab,0,4,table.scale,x0,&y0,&yp);
+ evaluate_table(table.data,0,4,table.scale,x0,&y0,&yp);
fprintf(fp,"%15.10e %15.10e %15.10e\n",x0,y0,yp);
}
Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
- evaluate_table(table.tab,0,4,table.scale,r,&y0,&yp);
+ evaluate_table(table.data,0,4,table.scale,r,&y0,&yp);
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);
abs_error_r=fabs(y0-Vtab);
* to do this :-)
*/
- snew_aligned(table.tab,4*nx,16);
-
- copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,table.tab);
+ snew_aligned(table.data,4*nx,16);
+
+ copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,1.0,table.data);
if(bDebugMode())
{
{
/* x0=i*table.r/(5*table.n); */
x0 = i*table.r/(5*(table.n-1));
- evaluate_table(table.tab,0,4,table.scale,x0,&y0,&yp);
+ evaluate_table(table.data,0,4,table.scale,x0,&y0,&yp);
fprintf(fp,"%15.10e %15.10e %15.10e\n",x0,y0,yp);
}
done_tabledata(&(td[0]));
sfree(td);
+
+ table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
+ table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
+ table.formatsize = 4;
+ table.ninteractions = 3;
+ table.stride = table.formatsize*table.ninteractions;
+
return table;
}
}
tab.n = td.nx;
tab.scale = td.tabscale;
- snew(tab.tab,tab.n*4);
- copy2table(tab.n,0,4,td.x,td.v,td.f,tab.tab);
+ snew(tab.data,tab.n*4);
+ copy2table(tab.n,0,4,td.x,td.v,td.f,1.0,tab.data);
done_tabledata(&td);
return tab;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff