A problem with lj/smooth/linear potential

Hi there,
I used lj/smooth/linear potential in my MD simulation by Lammps, but I found that the total potential energy of my system is different when different number of processors is specified. The input script and thermo output was attached in my previous post just several posts prior to this one: “A couple of problems about compute group/group command”. I shall attach them again if you could not see the prior posts. May I ask if there are some suggestions? Thank you very much!

steffani,

Hi there,
     I used lj/smooth/linear potential in my MD simulation by Lammps, but I
found that the total potential energy of my system is different when
different number of processors is specified. The input script and thermo
output was attached in my previous post just several posts prior to this
one: "A couple of problems about compute group/group command". I shall
attach them again if you could not see the prior posts. May I ask if there
are some suggestions? Thank you very much!

first of all. the input you posted is *far* to complex
to conclude anything from it.

if you want anybody to take a serious look at the issue,
you have to boil it down to a "minimal" input with just
a small number of atoms, and perhaps only two or
even one layer.

i also would strongly recommend to try without
pair_style hybrid. i have a strange suspicion, that
there is something else going on that is very different
from your speculations.

cheers,
     axel.

Hi Axel,
thank you for your suggestions about tricks of posting. I reduced the system to 6 layers and the input script is reduced as follows, which I believe is simple enough for a serious look:

group layer3 type <> 3 3
group layer4 type <> 4 4
group layer1 type <> 1 1
group layer6 type <> 6 6

pair_style hybrid lj/smooth/linear 7.13176 tersoff
pair_coeff * * tersoff SiC.tersoff C C C C C C
pair_coeff 1 26 lj/smooth/linear .003 3.4
pair_coeff 2 3
6 lj/smooth/linear .003 3.4
pair_coeff 3 46 lj/smooth/linear .003 3.4
pair_coeff 4 5
6 lj/smooth/linear .003 3.4
pair_coeff 5 6 lj/smooth/linear .003 3.4

neighbor 2.0 bin
neigh_modify delay 0 every 20 check no

fix 1 all npt temp 300.0 300.0 0.1 x 1.0 1.0 0.2 y 1.0 1.0 0.2 z 1.0 1.0 0.2 drag 0.2
compute PE34 layer3 group/group layer4
compute PE16 layer1 group/group layer6
fix calcPELJ all ave/time 1 500 500 c_PE34 c_PE16 file PELJ.profile
thermo_style custom step temp pe ke etotal c_PE34 c_PE16
thermo 100

run 10000

Here is the thermo output for 1 processor run:

Step Temp PotEng KinEng TotEng PE34 PE16
0 300 -12359.951 65.108335 -12294.842 -8.6168987 -20.173609
100 158.22946 -12457.059 34.340189 -12422.719 -7.8785547 -17.444566
200 157.53566 -12453.449 34.189616 -12419.259 -7.8959591 -17.447109

Step Temp PotEng KinEng TotEng PE34 PE16
0 300 -12356.539 65.108335 -12291.431 -8.6168987 -20.173609
100 158.22946 -12456.527 34.340189 -12422.187 -7.8785547 -17.444566
200 157.53566 -12453.018 34.189616 -12418.829 -7.8959591 -17.447109

It is still shown that the potential energy of the system is different for different number of processors.
This is the smallest structure to boil down to based on the current cutoff. And I cannot avoid using hybrid because I have to hold the atoms in plane by tersoff potential. Besides, based on the fact that every result remains exactly the same no matter how many processors are specified when LJ/smooth/linear is replaced by morse potential, I don’t see if there is anything wrong with hybrid.

Sorry, I forgot to say that the second thermo output is for the 12 processors run.

Hi Axel,
thank you for your suggestions about tricks of posting. I reduced the system to 6 layers and the input script is reduced as follows, which I believe is simple enough for a serious look:

group layer3 type <> 3 3
group layer4 type <> 4 4
group layer1 type <> 1 1
group layer6 type <> 6 6

pair_style hybrid lj/smooth/linear 7.13176 tersoff
pair_coeff * * tersoff SiC.tersoff C C C C C C
pair_coeff 1 26 lj/smooth/linear .003 3.4
pair_coeff 2 3
6 lj/smooth/linear .003 3.4
pair_coeff 3 46 lj/smooth/linear .003 3.4
pair_coeff 4 5
6 lj/smooth/linear .003 3.4
pair_coeff 5 6 lj/smooth/linear .003 3.4

neighbor 2.0 bin
neigh_modify delay 0 every 20 check no

fix 1 all npt temp 300.0 300.0 0.1 x 1.0 1.0 0.2 y 1.0 1.0 0.2 z 1.0 1.0 0.2 drag 0.2
compute PE34 layer3 group/group layer4
compute PE16 layer1 group/group layer6
fix calcPELJ all ave/time 1 500 500 c_PE34 c_PE16 file PELJ.profile
thermo_style custom step temp pe ke etotal c_PE34 c_PE16
thermo 100

run 10000

Here is the thermo output for 1 processor run:

Step Temp PotEng KinEng TotEng PE34 PE16
0 300 -12359.951 65.108335 -12294.842 -8.6168987 -20.173609
100 158.22946 -12457.059 34.340189 -12422.719 -7.8785547 -17.444566
200 157.53566 -12453.449 34.189616 -12419.259 -7.8959591 -17.447109

Here is the thermo output for 1 processor run:

Hi Axel,
      thank you for your suggestions about tricks of posting. I reduced the
system to 6 layers and the input script is reduced as follows, which I

why not 2?

believe is simple enough for a serious look:

from what you provided, nobody can run.
unless you provide a *complete* set of inputs
and reduce it to the smallest unit that reproduces
the issue, chances are small that anybody
will bite. the easier you make it for somebody
to reproduce your numbers, they more likely,
this person will help you.

the time it takes me, for example, to come up
with a way to construct your geometry without
really knowing what you have, and then strip
it down to the essentials, i could have answered
quite a few e-mails from other people.

axel.

Why not 2?
If your cutoff is 7.**, and if you use periodic boundary condition and you make your boxsize <7, it’ll cause a problem. You will find it in any MD book.

If you have to run simulations of everyone with a problem to tackle the problems, I can send you my data file together with the script.

Why not 2?
If your cutoff is 7.**, and if you use periodic boundary condition and you
make your boxsize <7, it'll cause a problem. You will find it in any MD
book.

Then you don't really understand how LAMMPS works. You can still run
a box less than 7 if your cutoff is 7 with LAMMPS.

If you have to run simulations of everyone with a problem to tackle the
problems, I can send you my data file together with the script.

You have to help people who are trying to help you, especially for
free and on their own time.

Please first see if you can reproduce the difference using different
number of cores with just pair_style lj/smooth/linear ( no hybrid with
Tersoff). Like Axel mentioned, problem may come from hybrid.

Thank you, Ray. Sure, I would be equally happy if I or someone else could have the problem solved. I appreciate you and Dr. Axel for your helps. And I would do anything I could to help as well. If I caused any offense due to misunderstanding, I apologize.

I saw before that lammps still runs when the cutoff is larger than half the box size. But I need to evaluate the interaction energy between layers, so I am not confident that too small a boxsize will not cause a problem. May I ask how lammps manage to run with box size smaller than twice cutoff? Thank you very much!

I am still conservative about the blame to hybrid, but I will do test runs for sure and I would only make comments based on facts.

Steffani,
For debugging it does not matter to have a meaningful input. It simply needs to reproduce the inconsistency. Don’t rule out using hybrid. A lot of problems have happened because people underestimated the impact of manybody terms that were included in ways they did not expect.
Axel.

Don't rule out using hybrid. A lot of
problems have happened because people underestimated the impact of manybody
terms that were included in ways they did not expect.

I want to emphasize this point. I don't understand your model, using hybrid,
defining Tersoff for all interactions, then defining some with
lj/smooth. B/c you are using hybrid this will have the effect of
overwriting some of the Tersoff interactions with the lj/smooth ones. That
is an inconsistent thing to do with a manybody potential, b/c of the way
it finds and computes interactions, e.g. 3 body terms. It's possible that
is just a wrong thing to do, and when the nieighbor lists for Tersoff
are formed,
the results differ depending on what order the atoms and their types are listed.

Hence different answers on different #s of processors.

Steve