problem on using "pair_type table" and "pair_write"

Dear all,

I met one question when using “pair_type table” and “pair_write”, which stuck me for a long time:

the related part in my input file is this:
pair_type table linear 500
pair_coeff 1 1 pair11.table ENTRY
pair_write 1 1 500 r 0.0001 1.9671 table.txt Yukawa

the content of pair11.table is:
#pair table
ENTRY
N 500 RSQ 0.0001 1.9671

1 1E-8 24.83089 -89.4697
2 1.63368E-5 24.47838 -89.4697
3 6.37406E-5 24.12587 -89.3575
4 1.42222E-4 23.77335 -89.3575
5 2.51778E-4 23.42084 -89.4694
6 3.92412E-4 23.06833 -89.4694
7 5.64124E-4 22.71582 -89.4709
8 7.66913E-4 22.3633 -89.3575
9 0.001 22.01079 -89.3563
10 0.00127 21.65828 -89.4694
11 0.00156 21.30577 -89.4706
.
.
.
498 3.83854 -0.00195 -0.0063
499 3.854 -0.00198 -0.0063
500 3.86948 -0.002 -0.0051

According to the lammps mannual, "if you want the interpolation tables of length Ntable to match exactly what is in the tabulated file (with effectively no preliminary interpolation), you should set Ntable = Nfile, and use the “RSQ” or “BITMAP” parameter ", I use RSQ here. RSQ is followed by 2 values rlo (0.0001) and rhi(1.9671). The “rlo” should be 0.0 in my simulation (DPD simulation), but 0.0 is not allowed, so I choose 0.0001 instead of 0.0. The first column is the index from 1 to 500, the second column is ranged from rlorlo (1E-8) to rhirhi(3.86948), the third column is the energy and the forth column is the force which is obtained by using central difference method on energy and distance.

But after using “pair_write”, I find the content in table.txt is:

1 1E-4 27.3304 -84.4702
2 0.00404 27.31 -3603.68
3 0.00798 27.2882 -7079.45
4 0.01193 27.2652 -10469.4
5 0.01587 27.2407 -13731.3
6 0.01981 27.215 -16822.7
7 0.02375 27.1879 -19701.2
8 0.02769 27.1595 -22324.5
9 0.03164 27.1297 -24650.3
10 0.03558 27.0986 -26636.2
.
.
.
498 1.95922 -0.0019 -0.0051
499 1.96316 -0.0019 -0.0063
500 1.9671 0 0

Comparing it to the input table file(pair11.table, discussed aboved ), the third column is the energy starting from 27.3304 which is different from that in pair11.table starting from 24.83089, the fourth colume is the force which show the worst comparison because the force increases to very huge values in the short range.

To correct this, I tried many ways suggest by the lammps manual:

  • Vary the number of table points; you may need to use more than you think to get good resolution.
  • Always use the pair_write command to produce a plot of what the final interpolated potential looks like. This can show up interpolation “features” you may not like.
  • Start with the linear style; it’s the style least likely to have problems.
  • Use N in the pair_style command equal to the “N” in the tabulation file, and use the “RSQ” or “BITMAP” parameter, so additional interpolation is not needed. See discussion below.
  • Use as large an inner cutoff as possible. This avoids fitting splines to very steep parts of the potential.

I tried suggest 1, I changed the table points from N=500 to N=10000, but it did not work for my case. For suggest 3, I indeed used the linear style. For suggest 4, I indeed used the N in the pair_style command equal to the “N” in the tabulation file which is 500, and use the “RSQ” parameter. For suggest 5, I cannot use a too large inner cutoff in DPD simulation, because the particle in DPD simulation can even overlap (r=0.0). If I use a large inner cutoff, the error"Pair distance < table inner cutoff" will show up after a few time steps.

I also looked into the codes of pair_table.cpp and pair_table.h in directory of src in lammps to see if I can find what is wrong with my problem. I noticed that when “match=1 if don’t need to spline read-in tables” and the corresponding codes are:
if (tabstyle==LINEAR && tb->ninput==tablength && tb->rflag==RSQ && rhi==tb->cut) tb->match=1

My case satisfies all the condition it lists and there should be match=1, therefore the read-in table (pair11.table) should not be splined. But the results are beyond my expection. I have no idea why the energy and force (especially the force) in output table file is so different from the read-in table at the short distance. I would like to know how I can get the same output table file as that of input table file in my case. Any help and suggestions will be greatly appriciated.

Regards,
Yingdong

I'll see if Paul can look into this.
He made one recent bug fix with respect
to the slope at end points for the splining.
Putting a few print statements in the code
should verify whether splining is actually
happening or not, since you are trying
to avoid it altogether.

Steve

Yes, I'd like to take a closer look at this.

Yingdong: could you send me your "pair11.table" file?

Thanks,

Paul

Dear Paul,

The input “pair11.table” together with the output “table.txt” are attached in this e-mail. Thanks so much! Have a great day!

Regards,
Yingdong

pair11.table (13.5 KB)

table.txt (15.2 KB)

Your table entries are evenly spaced in R instead of RSQ, so LAMMPS can’t use your values as-is and has to spline.

See: http://lammps.sandia.gov/doc/pair_table.html

Specifically:

“If used, the parameters “R” or “RSQ” are followed by 2 values rlo and rhi. If specified, the distance associated with each energy and force value is computed from these 2 values (at high accuracy), rather than using the (low-accuracy) value listed in each line of the table. For “R”, distances uniformly spaced between rlo and rhi are computed; for “RSQ”, squared distances uniformly spaced between rlo*rlo and rhi*rhi are computed.”

So I’d recommend spacing them evenly in RSQ and trying again.

Paul