# How to define Potential in tabular form

Hi all,

In a paper, the author gave some graph of energy versus angle, distance etc and those graphs have no single peak distribution. They mentioned that they defined the potential in LAMMPS using tabular form as there is no single peak. Would you please help me to know how can I define potential in LAMMPS in tabular form using those graphical values.

Thanks
Best Regards

Salah

hello,

you should try to read the manual first. if you do, you’d realize you should try the pair_style table and construct
the table accordingly http://lammps.sandia.gov/doc/pair_table.html

Remember to always check afterwards with pair_write command, so you can see if the table is correct

pablo

What is “single peak distribution”?

A potential, whether in tabular or analytical form, if correctly tabulated should not have anything to do with its potential energy surface.

pair_style table explains how to tabulate a potential.

Ray

Hi all,

In a paper, the author gave some graph of energy versus angle, distance etc and those graphs have no single peak distribution. They mentioned that they defined the potential in LAMMPS using tabular form as there is no single peak. Would you please help me to know how can I define potential in LAMMPS in tabular form using those graphical values.

do you want to implement the exact same potential of the publication
you are referring to, or just any tabulated potential?

for the second case, you already have some suggestions, for the first
case, i would contact the authors and ask for copies, or check, if
those tables. if that doesn't work, you can always get back to plan b
and tabulate yourself...

axel.

Hi,

I have defined the potential table for angle potential successfully but in case of dihedral it is showing " Invalid dihedral style (.. /force.cpp:406)" and from the LAMMPS manual the reason is mentioned " The choice of dihedral style is unknown"

I am using LAMMPS Windows serial executable --- C++ version (12 April 2013)

I have used the following command:
dihedral_style table spline 400
dihedral_coeff 1 pottorsion.txt torsional

Thanks
Best Regards

Md Salah Uddin

Hi,

I have defined the potential table for angle potential successfully but in case of dihedral it is showing " Invalid dihedral style (.. /force.cpp:406)" and from the LAMMPS manual the reason is mentioned " The choice of dihedral style is unknown"

that is correct.

I am using LAMMPS Windows serial executable --- C++ version (12 April 2013)
I have used the following command:
dihedral_style table spline 400
dihedral_coeff 1 pottorsion.txt torsional

please look carefully at the documentation of dihedral_style table and
the description of the windows executable.
the style is part of the USER-MISC package, but the windows executable
that you have doesn't include this package.

you can try using a windows executable from here

http://git.icms.temple.edu/rpm/windows.html

which *does* include USER-MISC and thus this style (and its

axel.

Hi,
Thanks a lot. That problem has been solved and encountering another problem regarding pair coefficient. I am trying to solve it myself first. If I can't, I will ask you again. Appreciate your help.

Thanks
Best Regards

Md Salah Uddin

Hi,

This time it is showing illegal pair coefficient command. I have checked the syntax numerous time. I have tried with different options like spline, linear with optional parameter R/RSQ again with same N for the pair style command and tabulation file. But nothing is working. Then I tried to add some cutoff distance with different values from 1 to 11 but every time it is showing numeric index is out of bound. My input code and the table file are as follows (the resolution is very high but due should this error come):

pair_style table spline 400
pair_coeff 1 potpair.txt couple

#pair potential
couple
N 44 R 2.4 11.0

1 2.4 3.066020 2.09755
2 2.6 2.646510 2.47535
3 2.8 2.151440 1.2167
4 3 1.908100 1.1929
5 3.2 1.669520 0.72005
6 3.4 1.525510 0.3245
7 3.6 1.460610 0.5601
8 3.8 1.348590 0.20605
9 4 1.307380 0.17055
10 4.2 1.273270 0.15875
11 4.4 1.241520 0.1528
12 4.6 1.210960 0.1587
13 4.8 1.179220 0.1824
14 5 1.142740 0.20015
15 5.2 1.102710 0.20605
16 5.4 1.061500 0.2238
17 5.6 1.016740 0.211965
18 5.8 0.974347 0.158745
19 6 0.942598 0.175535
20 6.2 0.907491 0.07593
21 6.4 0.892305 0.046365
22 6.6 0.883032 0.0109
23 6.8 0.880852 0.022725
24 7 0.876307 0.04639
25 7.2 0.867029 0.05229
26 7.4 0.856571 0.052285
27 7.6 0.846114 0.046375
28 7.8 0.836839 0.022725
29 8 0.832294 0.004995
30 8.2 0.831295 -0.00092
31 8.4 0.831479 0 .016815
32 8.6 0.828116 0.028645
33 8.8 0.822387 0.02865
34 9 0.816657 0.02272
35 9.2 0.812113 0.030565
36 9.4 0.806000 -0.00784
37 9.6 0.807568 0.016825
38 9.8 0.804203 0.010905
39 10 0.802022 0.022725
40 10.2 0.797477 0.022735
41 10.4 0.792930 0.0109
42 10.6 0.790750 -0.006825
43 10.8 0.792115 0.004995
44 11 0.791116 0.004995

Thanks
Best Regards

Md Salah Uddin

potpair.txt (1.11 KB)

Hi,

This time it is showing illegal pair coefficient command. I have checked the syntax numerous time. I have tried with different options like spline, linear with optional parameter R/RSQ again with same N for the pair style command and tabulation file. But nothing is working. Then I tried to add some cutoff distance with different values from 1 to 11 but every time it is showing numeric index is out of bound. My input code and the table file are as follows (the resolution is very high but due should this error come):

pair_style table spline 400
pair_coeff 1 potpair.txt couple

docs again. the single number only applies to data files.

axel.

Thanks. Now that is gone. But is now showing "pair distance <table inner cutoff>" that means two atoms are closer together than the pairwise table allows. Do I need to modify my model then or is there any other options?

Best Regards

Md Salah Uddin

Thanks. Now that is gone. But is now showing "pair distance <table inner cutoff>" that means two atoms are closer together than the pairwise table allows. Do I need to modify my model then or is there any other options?

what do *you* think? it is *your* simulation. it is *your* data. you
figure it out.

axel.

Good Morning Dr. Axel Kohlmeyer,

I have modified my polymer model and rearranged all the atoms which seems very closer. But still it is showing "pair distance <table inner cutoff>". Would you please give me some hints what is the minimum and maximum distance between two atoms that pair table allows? Please give me some idea how to solve this.

Thanks
Best Regards

Md Salah Uddin

Well, I am trying to do the simulation of a single chain to avoid the cut off distance problem temporarily. I have a chain of 25 atoms and i am pulling it in x direction using deform erate. But I visualize the results, initially the chain is divided into two and the two splits are parallel each other apparently. and I am not getting any pressure tensor values. The two separate portions are going away after some time steps. Would you please help me regarding this?

Thanks
Best Regards
Salah

Well, I am trying to do the simulation of a single chain to avoid the cut off distance problem temporarily. I have a chain of 25 atoms and i am pulling it in x direction using deform erate. But I visualize the results, initially the chain is divided into two and the two splits are parallel each other apparently. and I am not getting any pressure tensor values. The two separate portions are going away after some time steps. Would you please help me regarding this?

my guess is that you must be doing something wrong.

axel.

Would you mind if I send you my input and data files and request you to look at it? It is just about 25 atoms of a single polymer chain. I am new at MD simulation and I am really stuck at this point. If you have no problem, please let me know then I will give you my files.

Thanks
Best Regards

Md Salah Uddin

couple
N 44 R 2.4 11.0

1 2.4 3.066020 2.09755
2 2.6 2.646510 2.47535
3 2.8 2.151440 1.2167
4 3 1.908100 1.1929
5 3.2 1.669520 0.72005
6 3.4 1.525510 0.3245
7 3.6 1.460610 0.5601
8 3.8 1.348590 0.20605
9 4 1.307380 0.17055
10 4.2 1.273270 0.15875
11 4.4 1.241520 0.1528

It looks like you omitted the part of the table which describes the
repulsive core, so it is not surprising the atoms will come close
together. The pairwise-force table begins at 2.4, and the energy
between the atoms at this point is not very large. (details: Assuming
these energies are in kCal/mole (default units), then at room
temperature, these energies are accessable, -especially if you begin
the simulation with the atomic positions not well equilibrated.)

Why don't you create a new table, starting with a smaller radius?
(...and using a higher density of entries. It looks like the force
and energies change a great deal every 0.2 Angstroms. tables often
have a few hundred entries or more.)

Cheers
Andrew

P.S.
(unrelated comment), be sure to use the "pair_write" command to make
sure there are no numerical instabilities in your pairwise forces and
energies

Thanks a lot. I will try this way now. I appreciate your suggestions.

Best Regards
Md Salah Uddin

Good Morning,

Using Fix NVE/limit I can supersede that error "Pair table inner cut off distance" and I am using fix deform as follows:

fix 1 all nve/limit 0.05
fix 3 all langevin 298 298 1 290
fix 2 all deform 1 x erate 1e-5 units box remap x

Taking the pressure tensor as stress I am getting quite reasonable stress-strain curve. Using Fix NVE/limit, can I accept this stress-strain values?

Thanks
Best Regards

Md Salah Uddin

Good Morning,

Using Fix NVE/limit I can supersede that error "Pair table inner cut off distance" and I am using fix deform as follows:

fix 1 all nve/limit 0.05
fix 3 all langevin 298 298 1 290
fix 2 all deform 1 x erate 1e-5 units box remap x

Taking the pressure tensor as stress I am getting quite reasonable stress-strain curve. Using Fix NVE/limit, can I accept this stress-strain values?

Personally speaking, it's okay with me. I will leave it to you to
check that the langevin damping rates ("1") and deformation rates
(erate 1e-5) are reasonable.

potential energy or free energy?

I would like to point out that you are running these simulations at
non-zero temperature (298K). This means that the potential energy
will fluctuate randomly. It also means that the free_energy !=
potential_energy. So plotting potential energy as a function of
strain is not really correct. Formally speaking, you should plot
free_energy as a function of strain. Unfortunately calculating
free_energy is more difficult. (In my field, soft biomaterials, there
are many ways to do this: free energy perturbation, thermodynamics
integration, umbrella-sampling & WHAM, and other non-equilibrium
methods based on the Jarzynski equality/Crooke's Theorem.)

However, hopefully you can avoid that extra work. A temperature of
298K is low, and perhaps the stress-strain curve for your material
does not depend on temperature much at these low (298K) temperatures.
In that case, run Langevin dynamics at zero temperature (or a
temperature close to zero). For example:

fix 3 all langevin 0.001 0.001 1 290

That way you can plot potential energy as a function of strain.
I am blindly guessing for metals and semiconductors (for example) the
difference between free energy and energy (potential+kinetic) is
negligible at these low temperatures. Please check with someone who
is in your field of research who is more knowledgeable than me. How
do they do it? Good luck.

Andrew