Project

General

Profile

gromacs-4.5.3-gbsalt.patch

Patch for gbsalt term against GROMACS 4.5.3 - Matthew Zwier, 02/26/2013 05:33 PM

View differences:

gromacs-4.5.3-gbsalt/include/force.h 2011-03-02 09:30:56.000000000 -0500
114 114
                                  const char *fn,
115 115
                                  real rtab);
116 116

  
117
t_forcetable make_gb_table2(FILE *out,const output_env_t oenv,
118
                                  const t_forcerec *fr,
119
                                  const char *fn,
120
                                  real rtab);
121

  
117 122
void pr_forcerec(FILE *fplog,t_forcerec *fr,t_commrec *cr);
118 123

  
119 124
void
gromacs-4.5.3-gbsalt/include/genborn.h 2011-03-02 09:30:56.000000000 -0500
71 71
				  real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
72 72
				  real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph);
73 73

  
74

  
74
real gb_bonds_tab_salt(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
75
                  real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
76
                  real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph, 
77
                  real salt_conc, real *GBtab2, real *p_gbtabscale2);
75 78

  
76 79
void gb_pd_send(t_commrec *cr, real *send_data, int nr);
77 80

  
gromacs-4.5.3-gbsalt/include/types/forcerec.h 2011-03-02 09:30:56.000000000 -0500
299 299
  real gb_epsilon_solvent;
300 300
  /* Table data for GB */
301 301
  t_forcetable gbtab;
302
  /* Alternate table to include salt concentration term */
303
  t_forcetable gbtab2; 
302 304
  /* VdW radius for each atomtype (dim is thus ntype) */
303 305
  real *atype_radius;
304 306
  /* Effective radius (derived from effective volume) for each type */
......
316 318
  real gbtabscale;
317 319
  /* Table range for GB */
318 320
  real gbtabr;
321

  
322
  /* Table scale for GB */
323
  real gbtabscale2;
324
  /* Table range for GB */
325
  real gbtabr2;
326

  
319 327
  /* GB neighborlists (the sr list will contain for each atom all other atoms                                            
320 328
   * (for use in the SA calculation) and the lr list will contain                                                  
321 329
   * for each atom all atoms 1-4 or greater (for use in the GB calculation)                                        
gromacs-4.5.3-gbsalt/src/mdlib/forcerec.c 2011-03-02 09:30:57.000000000 -0500
1598 1598
	{
1599 1599
#ifdef GMX_DOUBLE
1600 1600
		fr->gbtabscale=2000;
1601
        fr->gbtabscale2=20000;
1601 1602
#else
1602 1603
		fr->gbtabscale=500;
1604
        fr->gbtabscale2=5000;
1603 1605
#endif
1604 1606
		
1605
		fr->gbtabr=100;
1607
		fr->gbtabr=100;        
1606 1608
		fr->gbtab=make_gb_table(fp,oenv,fr,tabpfn,fr->gbtabscale);
1609
        fr->gbtabr2=10;
1610
        fr->gbtab2=make_gb_table2(fp,oenv,fr,tabpfn,fr->gbtabscale2);
1607 1611

  
1608 1612
        init_gb(&fr->born,cr,fr,ir,mtop,ir->rgbradii,ir->gb_algorithm);
1609 1613

  
gromacs-4.5.3-gbsalt/src/mdlib/genborn.c 2011-03-02 09:30:57.000000000 -0500
1330 1330
    return 0;        
1331 1331
}
1332 1332

  
1333
real gb_bonds_tab_salt(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1334
                  real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1335
                  real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph, 
1336
                  real salt_conc, real *GBtab2, real *p_gbtabscale2)
1337
{
1338
    int i,j,n0,m,nnn,type,ai,aj;
1339
	int ki;
1340
 
1341
	real isai,isaj;
1342
    real r,rsq11;
1343
    real rinv11,iq;
1344
    real isaprod,gbscale,gbtabscale,Y,F,Geps,Heps2,Fp,VV,rt,eps,eps2,FF;
1345
    //real qq,fijC;
1346
    real vgb,fgb,vcoul,dvdatmp,fscal,dvdaj;
1347
    real vctot;
1333 1348

  
1349
    real invxi, diffinvxi, u, saltterm, expterm, ut, ueps, ueps2, invu, diffu, kqq;
1350
    real gbtabscale2;
1351
    int un0, unnn;
1352
    real kappa, inv_gb_epsilon_solvent, inv_epsilon_r;
1353
	rvec dx;
1354
	ivec dt;
1355
	
1356
    kappa = 0.316*sqrt(salt_conc);
1357
    inv_gb_epsilon_solvent = 1.0/gb_epsilon_solvent;
1358
    inv_epsilon_r = 1.0/epsilon_r;
1359
    t_iatom *forceatoms;
1360

  
1361
    gbtabscale=*p_gbtabscale;
1362
    gbtabscale2=*p_gbtabscale2;
1363
    vctot = 0.0;
1364
    
1365
    for(j=F_GB12;j<=F_GB14;j++)
1366
    {
1367
        forceatoms = idef->il[j].iatoms;
1368
        
1369
        for(i=0;i<idef->il[j].nr; )
1370
        {
1371
            /* To avoid reading in the interaction type, we just increment i to pass over
1372
             * the types in the forceatoms array, this saves some memory accesses
1373
             */
1374
            i++;
1375
            ai            = forceatoms[i++];
1376
            aj            = forceatoms[i++];
1377
			
1378
			ki            = pbc_rvec_sub(pbc,x[ai],x[aj],dx);           
1379
			rsq11         = iprod(dx,dx); //rij^2
1380
			
1381
			isai          = invsqrta[ai]; // 1/sqrt(Ri)
1382
			iq            = (-1)*facel*charge[ai];
1383
			
1384
            rinv11        = gmx_invsqrt(rsq11); // 1/rij
1385
            isaj          = invsqrta[aj]; // 1/sqrt(Rj)
1386
            isaprod       = isai*isaj;    // 1/sqrt(Ri*Rj)
1387
            //qq            = isaprod*iq*charge[aj]; 
1388
            gbscale       = isaprod*gbtabscale; //500/sqrt(Ri*Rj) single pres
1389
            r             = rsq11*rinv11; //rij
1390
            rt            = r*gbscale;    //location of distance in GBTab
1391
            n0            = rt;           //implicit cast to int
1392
            eps           = rt-n0;        //fractional part of rt (n0 + eps = rt) see 6.7 in the GMX manual
1393
            eps2          = eps*eps;      //u^2
1394
            nnn           = 4*n0;         
1395
            Y             = GBtab[nnn];
1396
            F             = GBtab[nnn+1];
1397
            Geps          = eps*GBtab[nnn+2];
1398
            Heps2         = eps2*GBtab[nnn+3];
1399
            Fp            = F+Geps+Heps2;
1400
            VV            = Y+eps*Fp;          //GBtab[nnn] + u*GBtab[nnn+1] + u^2*GBtab[nnn+2] + u^3*GBtab[nnn+3]
1401
            FF            = Fp+Geps+2.0*Heps2; //GBtab[nnn+1] + 2*u*GBtab[nnn+2] + 3u^2*GBtab[nnn+3]
1402
            //vgb           = qq*VV;             
1403
            //fijC          = qq*FF*gbscale;
1404

  
1405
            invxi         = 1.0 / VV; /* 1.0/xi(x) -> xi(x) = 1.0/sqrt(x^2 + exp(-0.25*x^2)) */
1406
            /* Now determine the e(-u)/u value for this atom pair */
1407
            u             = kappa * invxi / isaprod; /* kappa * sqrt(Ri*Rj) / xi(x) */
1408
            ut            = u*gbtabscale2;
1409
            un0           = ut;
1410
            ueps          = ut-un0;
1411
            ueps2         = ueps*ueps;
1412
            unnn          = 4*un0;      
1413
            expterm       = GBtab2[unnn] + ueps*GBtab2[unnn+1] + ueps2*GBtab2[unnn+2] + ueps*ueps2*GBtab2[unnn+3]; /* exp(-u)/u */
1414
            //printf("u:%.16f e^(-u)/u:%.16f kappa:%.16f 1/Ri:%.16f 1/Rj:%.16f xi:%.16f x:%.16f\n",u,expterm,kappa,isai,isaj,VV,r*isaprod);
1415
            saltterm      = inv_gb_epsilon_solvent*expterm; /* (1.0/epsilon_w)*(exp(-u)/u) */
1416
            invu          = 1.0/u;
1417
            kqq           = kappa*iq*charge[aj]; /* kappa*qi*qj  */
1418
            diffu         = -inv_epsilon_r*invu*invu+saltterm*(invu+1.0); /* dvgb/du */
1419
            diffinvxi     = -invxi*invxi*FF*gbtabscale; /* d(invxi)/dx = 1.0/(xi(x)**2.0) * dxi/dx */ 
1420
            //printf("dinvxi:%.16f %.16f\n",diffinvxi,-FF*gbtabscale);
1421
            dvdatmp       = 0.5*kappa*diffu*kqq;
1422
            dvda[ai]      = dvda[ai] + dvdatmp*((isai/isaj)*invxi-diffinvxi*r*isai*isai); /* 0.5*kappa**2.0*qi*qj*dv/du*(sqrt(Rj)/(sqrt(Ri)*xi(x)) - d(invxi)/dx*rij/Ri) */
1423
            dvda[aj]      = dvda[aj] + dvdatmp*((isaj/isai)*invxi-diffinvxi*r*isaj*isaj); /* 0.5*kappa**2.0*qi*qj*dv/du*(sqrt(Ri)/(sqrt(Rj)*xi(x)) - d(invxi)/dx*rij/Rj) */
1424
            vgb           = kqq*(inv_epsilon_r*invu-saltterm); /* kappa*qi*qj*(1/(u*epsilon_r) - exp(-u)/(u*gb_epsilon_solvent)) */ 
1425
            vctot         = vctot + vgb;
1426
            fgb           = -kappa*kqq*diffu*diffinvxi*rinv11; /* -kappa**2.0*qi*qj*dv/du*(d(1/xi(x))/dx)/rij */
1427
			
1428
			if (graph) {
1429
				ivec_sub(SHIFT_IVEC(graph,ai),SHIFT_IVEC(graph,aj),dt);
1430
				ki=IVEC2IS(dt);
1431
			}
1432
			
1433
			for (m=0; (m<DIM); m++) {			/*  15		*/
1434
				fscal=fgb*dx[m];
1435
				f[ai][m]+=fscal;
1436
				f[aj][m]-=fscal;
1437
				fshift[ki][m]+=fscal;
1438
				fshift[CENTRAL][m]-=fscal;
1439
			}
1440
        }
1441
    }
1442
    
1443
    return vctot;
1444
}
1334 1445

  
1335 1446
real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1336 1447
                  real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
......
1371 1482
            aj            = forceatoms[i++];
1372 1483
			
1373 1484
			ki            = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
1374
			rsq11         = iprod(dx,dx);
1485
			rsq11         = iprod(dx,dx); //rij^2
1375 1486
			
1376
			isai          = invsqrta[ai];
1487
			isai          = invsqrta[ai]; // 1/sqrt(Ri)
1377 1488
			iq            = (-1)*facel*charge[ai];
1378 1489
			
1379
            rinv11        = gmx_invsqrt(rsq11);
1380
            isaj          = invsqrta[aj];
1381
            isaprod       = isai*isaj;
1382
            qq            = isaprod*iq*charge[aj];
1383
            gbscale       = isaprod*gbtabscale;
1384
            r             = rsq11*rinv11;
1385
            rt            = r*gbscale;
1386
            n0            = rt;
1387
            eps           = rt-n0;
1388
            eps2          = eps*eps;
1389
            nnn           = 4*n0;
1490
            rinv11        = gmx_invsqrt(rsq11); // 1/rij
1491
            isaj          = invsqrta[aj]; // 1/sqrt(Rj)
1492
            isaprod       = isai*isaj;    // 1/sqrt(Ri*Rj)
1493
            qq            = isaprod*iq*charge[aj]; 
1494
            gbscale       = isaprod*gbtabscale; //500/sqrt(Ri*Rj) single pres
1495
            r             = rsq11*rinv11; //rij
1496
            rt            = r*gbscale;    //location of distance in GBTab -> r/sqrt(RiRj)
1497
            n0            = rt;           //implicit cast to int
1498
            eps           = rt-n0;        //fractional part of rt (n0 + eps = rt) see 6.7 in the GMX manual
1499
            eps2          = eps*eps;      //u^2
1500
            nnn           = 4*n0;         
1390 1501
            Y             = GBtab[nnn];
1391 1502
            F             = GBtab[nnn+1];
1392 1503
            Geps          = eps*GBtab[nnn+2];
1393 1504
            Heps2         = eps2*GBtab[nnn+3];
1394 1505
            Fp            = F+Geps+Heps2;
1395
            VV            = Y+eps*Fp;
1396
            FF            = Fp+Geps+2.0*Heps2;
1397
            vgb           = qq*VV;
1398
            fijC          = qq*FF*gbscale;
1506
            VV            = Y+eps*Fp;          //GBtab[nnn] + u*GBtab[nnn+1] + u^2*GBtab[nnn+2] + u^3*GBtab[nnn+3]
1507
            FF            = Fp+Geps+2.0*Heps2; //GBtab[nnn+1] + 2*u*GBtab[nnn+2] + 3u^2*GBtab[nnn+3]
1508
            vgb           = qq*VV;             
1509
            fijC          = qq*FF*gbscale;     //gbtabscale*-facel*qiqj/(Ri*Rj) *(-de/dx) = (gbtabscale * qq / (Ri*Rj)) * (-de/dx)
1399 1510
            dvdatmp       = -(vgb+fijC*r)*0.5;
1400 1511
            dvda[aj]      = dvda[aj] + dvdatmp*isaj*isaj;
1401 1512
            dvda[ai]      = dvda[ai] + dvdatmp*isai*isai;
......
1470 1581
    
1471 1582
}
1472 1583

  
1584
real calc_gb_selfcorrections_salt(t_commrec *cr, int natoms, 
1585
                 real *charge, gmx_genborn_t *born, real *dvda, t_mdatoms *md, double facel,
1586
                 real salt_conc, real *GBtab2, real *p_gbtabscale2)
1587
{    
1588
    int i,ai,at0,at1;
1589
    real rai,e,derb,q,q2,fi,rai_inv,vtot;
1590
    real u, ut, ueps, ueps2, expterm, saltterm, kqq, inv_eru;
1591
    real gbtabscale2;
1592
    int un0, unnn;
1593
    real kappa, inv_gb_epsilon_solvent, inv_epsilon_r;	
1594

  
1595
    kappa = 0.316*sqrt(salt_conc);
1596
    inv_gb_epsilon_solvent = 1.0/born->gb_epsilon_solvent;
1597
    inv_epsilon_r = 1.0/born->epsilon_r;
1598

  
1599
    gbtabscale2 = *p_gbtabscale2;
1600

  
1601
    if(PARTDECOMP(cr))
1602
    {
1603
        pd_at_range(cr,&at0,&at1);
1604
    }
1605
    else if(DOMAINDECOMP(cr))
1606
    {
1607
        at0=0;
1608
        at1=cr->dd->nat_home;
1609
    }
1610
    else
1611
    {
1612
        at0=0;
1613
        at1=natoms;
1614
        
1615
    }
1616
        
1617
    vtot=0.0;
1618
    
1619
    /* Apply self corrections */
1620
    for(i=at0;i<at1;i++)
1621
    {
1622
        ai       = i;
1623
        
1624
        if(born->use[ai]==1)
1625
        {
1626
			rai           = born->bRad[ai];
1627
            rai_inv       = 1.0/rai;
1628
            q             = charge[ai];
1629
            q2            = q*q;
1630
            fi            = facel*q2;
1631
            kqq           = fi*kappa;
1632
            /* Now determine the e(-u)/u value for this atom pair */
1633
            u             = kappa * rai; /* kappa * Ri */            
1634
            ut            = u*gbtabscale2;
1635
            un0           = ut;
1636
            ueps          = ut-un0;
1637
            ueps2         = ueps*ueps;
1638
            unnn          = 4*un0;           
1639
            expterm       = GBtab2[unnn] + ueps*GBtab2[unnn+1] + ueps2*GBtab2[unnn+2] + ueps*ueps2*GBtab2[unnn+3]; /* exp(-u)/u */
1640
            saltterm      = expterm*inv_gb_epsilon_solvent; /* (1.0/epsilon_w)*(exp(-u)/u) */
1641
            inv_eru       = inv_epsilon_r/u; /* 1/(epsilon_r*u) */
1642
            e             = kqq*(inv_eru - saltterm); /* kappa*qi*qi(1.0/(epsilon_r*u)-(1.0/epsilon_w)*(exp(-u)/u)) */
1643
            dvda[ai]      -= 0.5*kqq*(-inv_eru*rai_inv+saltterm*(rai_inv + kappa));
1644
            vtot          -= 0.5*e;
1645
        }
1646
    }
1647
    
1648
   return vtot;    
1649
    
1650
}
1651

  
1473 1652
real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr,int natoms,gmx_genborn_t *born, gmx_localtop_t *top, 
1474 1653
                      const t_atomtypes *atype, real *dvda,int gb_algorithm, t_mdatoms *md)
1475 1654
{
......
1677 1856
  }
1678 1857
  
1679 1858
  /* Calculate the bonded GB-interactions using either table or analytical formula */
1680
    enerd->term[F_GBPOL]       += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1681
                                     fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1682
    
1683
    /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1684
    enerd->term[F_GBPOL]       += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);         
1859
    if(fr->userint1 == 0) //no salt
1860
    {
1861
        enerd->term[F_GBPOL]       += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1862
                                         fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1863
        /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1864
        enerd->term[F_GBPOL]       += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);       
1865
    }
1866
    else
1867
    {
1868
        enerd->term[F_GBPOL]       += gb_bonds_tab_salt(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1869
                                         fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph,
1870
                                         fr->userreal1, fr->gbtab2.tab, &(fr->gbtabscale2));
1871
        /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1872
        enerd->term[F_GBPOL]       += calc_gb_selfcorrections_salt(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac,
1873
                                         fr->userreal1, fr->gbtab2.tab, &(fr->gbtabscale2));
1874
    }
1875
  
1685 1876

  
1686 1877
    /* If parallel, sum the derivative of the potential w.r.t the born radii */
1687 1878
    if(PARTDECOMP(cr))
gromacs-4.5.3-gbsalt/src/mdlib/tables.c 2011-03-02 16:53:23.000000000 -0500
964 964
  return table;
965 965
}
966 966

  
967
t_forcetable make_gb_table2(FILE *out,const output_env_t oenv,
968
                           const t_forcerec *fr,
969
                           const char *fn,
970
                           real rtab)
971
{
972
	const char *fns[3] = { "gbctab_2.xvg", "gbdtab_2.xvg", "gbrtab_2.xvg" };
973
	const char *fns14[3] = { "gbctab14_2.xvg", "gbdtab14_2.xvg", "gbrtab14_2.xvg" };
974
	FILE        *fp;
975
	t_tabledata *td;
976
	gmx_bool        bReadTab,bGenTab;
977
	real        x0,y0,yp;
978
	int         i,j,k,nx,nx0,tabsel[etiNR];
979
	void *      p_tmp;
980
	double      r,r2,Vtab,Ftab,expterm,invr;
981
	
982
	t_forcetable table;
983
	
984
	double abs_error_r, abs_error_r2;
985
	double rel_error_r, rel_error_r2;
986
	double rel_error_r_old=0, rel_error_r2_old=0;
987
	double x0_r_error, x0_r2_error;
988
	
989
	
990
	/* Only set a Coulomb table for GB */
991
	/* 
992
	 tabsel[0]=etabGB;
993
	 tabsel[1]=-1;
994
	 tabsel[2]=-1;
995
	*/
996
	
997
	/* Set the table dimensions for GB, not really necessary to
998
	 * use etiNR (since we only have one table, but ...) 
999
	 */
1000
	snew(td,1);
1001
	table.r         = fr->gbtabr2;
1002
	table.scale     = fr->gbtabscale2;
1003
	table.scale_exp = 0;
1004
	table.n         = table.scale*table.r;
1005
	nx0             = 0;
1006
	nx              = table.scale*table.r;
1007
	
1008
	/* Check whether we have to read or generate 
1009
	 * We will always generate a table, so remove the read code
1010
	 * (Compare with original make_table function
1011
	 */
1012
	bReadTab = FALSE;
1013
	bGenTab  = TRUE;
1014
	
1015
	/* Each table type (e.g. coul,lj6,lj12) requires four 
1016
	 * numbers per datapoint. For performance reasons we want
1017
	 * the table data to be aligned to 16-byte. This is accomplished
1018
	 * by allocating 16 bytes extra to a temporary pointer, and then
1019
	 * calculating an aligned pointer. This new pointer must not be
1020
	 * used in a free() call, but thankfully we're sloppy enough not
1021
	 * to do this :-)
1022
	 */
1023
	
1024
	/* 4 fp entries per table point, nx+1 points, and 16 bytes extra 
1025
           to align it. */
1026
	p_tmp = malloc(4*(nx+1)*sizeof(real)+16);
1027
	
1028
	/* align it - size_t has the same same as a pointer */
1029
	table.tab = (real *) (((size_t) p_tmp + 16) & (~((size_t) 15)));  
1030
	
1031
	init_table(out,nx,nx0,table.scale,&(td[0]),!bReadTab);
1032
	
1033
	/* Local implementation so we don't have to use the etabGB
1034
	 * enum above, which will cause problems later when
1035
	 * making the other tables (right now even though we are using
1036
	 * GB, the normal Coulomb tables will be created, but this
1037
	 * will cause a problem since fr->eeltype==etabGB which will not
1038
	 * be defined in fill_table and set_table_type
1039
	 */
1040
	td->v[0] = td->f[0] = 1.0/0.0;
1041
	for(i=nx0+1;i<nx;i++)
1042
    {
1043
		Vtab    = 0.0;
1044
		Ftab    = 0.0;
1045
		r       = td->x[i];
1046
                invr    = 1.0/r;
1047
                expterm = exp(-r);
1048
		Vtab = expterm*invr;
1049
		Ftab = Vtab*(1.0+invr);
1050
		
1051
		/* Convert to single precision when we store to mem */
1052
		td->v[i]  = Vtab;
1053
		td->f[i]  = Ftab;
1054
		
1055
    }
1056
	
1057
	copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,table.tab);
1058
	
1059
	if(bDebugMode())
1060
    {
1061
		fp=xvgropen(fns[0],fns[0],"r","V",oenv);
1062
		/* plot the output 5 times denser than the table data */
1063
		/* for(i=5*nx0;i<5*table.n;i++) */
1064
		for(i=nx0;i<table.n;i++)
1065
		{
1066
			/* x0=i*table.r/(5*table.n); */
1067
			x0=i*table.r/table.n;
1068
			evaluate_table(table.tab,0,4,table.scale,x0,&y0,&yp);
1069
			fprintf(fp,"%15.10e  %15.10e  %15.10e\n",x0,y0,yp);
1070
			
1071
		}
1072
		ffclose(fp);
1073
    }
1074
	
1075
	/*
1076
	 for(i=100*nx0;i<99.81*table.n;i++)
1077
	 {
1078
	 r = i*table.r/(100*table.n);
1079
	 r2      = r*r;
1080
	 expterm = exp(-0.25*r2);
1081
	 
1082
	 Vtab = 1/sqrt(r2+expterm);
1083
	 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1084
	 
1085
	 
1086
	 evaluate_table(table.tab,0,4,table.scale,r,&y0,&yp);
1087
	 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);
1088
	 
1089
	 abs_error_r=fabs(y0-Vtab);
1090
	 abs_error_r2=fabs(yp-(-1)*Ftab);
1091
	 
1092
	 rel_error_r=abs_error_r/y0;
1093
	 rel_error_r2=fabs(abs_error_r2/yp);
1094
	 
1095
	 
1096
	 if(rel_error_r>rel_error_r_old)
1097
	 {
1098
	 rel_error_r_old=rel_error_r;
1099
	 x0_r_error=x0;
1100
	 }
1101
	 
1102
	 if(rel_error_r2>rel_error_r2_old)
1103
	 {
1104
	 rel_error_r2_old=rel_error_r2;
1105
	 x0_r2_error=x0;	
1106
	 }
1107
	 }
1108
	 
1109
	 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);
1110
	 printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1111
	 
1112
	 exit(1); */
1113
	done_tabledata(&(td[0]));
1114
	sfree(td);
1115
	
1116
	return table;
1117
	
1118
	
1119
}
1120

  
967 1121
t_forcetable make_gb_table(FILE *out,const output_env_t oenv,
968 1122
                           const t_forcerec *fr,
969 1123
                           const char *fn,
......
1048 1202
		
1049 1203
		Vtab = 1/sqrt(r2+expterm);
1050 1204
		Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1051
		
1205

  
1052 1206
		/* Convert to single precision when we store to mem */
1053 1207
		td->v[i]  = Vtab;
1054 1208
		td->f[i]  = Ftab;