pairwise interactions

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps users,
I wonder if any one can make some comments about these questions:
In calculating the pairwise interactions (in one of related cpp files such as pair_lj_cut_coul_long.cpp) why there is no periodic boundary condition applied in wrapping the distances? (I mean once the delx, dely, and delz btw i&j are calculated there is not any conditional statement such as if delx<-boxl/2 … & if delx>boxl/2…)

I have looked at all the pairs that pairwise interactions are calculated for in an example of a system including 300 atoms. I printed all the pairs that are listed (from the neighbor list) for the calculations,
Firstly, I found out there are atoms with ids higher than 299 in the list. (For instance atom number 500 or 600 or such). can anyone please tell me what are these?

Secondly, If I don’t use the neighbor list and simply loop over all atoms that are in cutoff (in a self-written program), let’s say I am finding all those that are in the attraction with atom No. “0”, there would many many more atoms!! why is this happening. If the neighbor list is making the code more efficient, it should not affect the pairs existed in the cutoff, right? otherwise the results will be different (as I found the forces different with lammps procedure and with my own simple example of looping over atoms as below)

for (int i=0; i<natom-1; i++)
for (int j=i; j<natom; j++)
delx=… dely=… delz…
if (delx> box/2)… if (delx<-box/2)… etc
rsq=delx^2+dely^2+delz^2
if (rsq< cutoffSqr)

calculations …

I appreciate any help/comments on these
emily

2

dear emily,

Dear Lammps users,
I wonder if any one can make some comments about these questions:
In calculating the pairwise interactions (in one of related cpp files such
as pair_lj_cut_coul_long.cpp) why there is no periodic boundary condition
applied in wrapping the distances? (I mean once the delx, dely, and delz btw

it is not needed. lammps uses a particle decomposition
and each node has "owned" and and "ghost" atoms. in
the case of just having a single domain (one processor),
those ghost atoms are the periodic images. this way
lammps is not bound in any way by minimum image
conventions and can handle cutoffs larger than half the
shortest box diameter.

i&j are calculated there is not any conditional statement such as *if
delx<-boxl/2* .... & *if delx>boxl/2...*)

I have looked at all the pairs that pairwise interactions are calculated for
in an example of a system including 300 atoms. I printed all the pairs that
are listed (from the neighbor list) for the calculations,
Firstly, I found out there are atoms with ids higher than 299 in the list.

yes, of course. those are the ghosts. see the variables atom->nlocal
and atom->nghost.

(For instance atom number 500 or 600 or such). can anyone please tell me
what are these?

see above.

i have to remark that it is very dangerous to look
at a code and not look at the documentation and
the papers that explain how it works.

Secondly, If I don't use the neighbor list and simply loop over all atoms
that are in cutoff (in a self-written program), let's say I am finding all
those that are in the attraction with atom No. "0", there would many many
more atoms!! why is this happening. If the neighbor list is making the code
more efficient, it should not affect the pairs existed in the cutoff, right?

all are there. you are not accounting for newton's
third law in your code apparently (and thus doing
twice as many force calculations).

mind you the "newton on/off" setting in LAMMPS applies to ghost
atoms only. for "local" atoms, the distinction is between a "half"
and a "full" neighbor list.

otherwise the results will be different (as I found the forces different
with lammps procedure and with my own simple example of looping over atoms
as below)

for (int i=0; i<natom-1; i++)
for (int j=i; j<natom; j++)
delx=... dely=... delz...
if (delx> box/2)... if (delx<-box/2)... etc
rsq=delx^2+dely^2+delz^2
if (rsq< cutoffSqr)

calculations ....

I appreciate any help/comments on these

ouch. this code fragment shows that you
still have to learn a lot about (efficient, parallel) MD.

you can find a quick summary of how to go from
a naiive MD implementation (like yours) to something
more efficient and parallel (up to using a cell list) in
this presentation (and also code examples for everything).

http://docs.google.com/viewer?a=v&pid=sites&srcid=ZGVmYXVsdGRvbWFpbnxha29obG1leXxneDo1OTQ0NTkzM2M1ZmY4ODA3&pli=1

and then i *strongly* suggest you read up the rest
about verlet lists and all the tricks of the trade in
a proper text book and steve's papers on the fundamental
ideas in the LAMMPS parallelization. failing to do
so will just make you misunderstand the code over
and over again.

cheers,
     axel.

3

Dear Axel,
Thanks for your detail comments and links, they were really useful to me.
I studied the ghost atom concept and the veret list and I am going to apply them. However, I want to make sure that the principles I am using work 100% correctly before making my code smarter.
Now, for a single processor the ghost atoms work the same as the periodic images. That means if I wrap the molecules dimension I must get those of ghost atoms (please correct me if I am wrong). I want to output the forces for 300 atoms (including all atoms and their equivalent ghosts). I mean, if I print (let’s say) f[0][0], it’s value show all the force counted for atom 0 either considered as owned or ghost with other ids (just want to see what the overall forces turned out to be)

Thanks again for you help
emily

4

Dear Axel,
Thanks for your detail comments and links, they were really useful to me.
I studied the ghost atom concept and the veret list and I am going to apply
them. However, I want to make sure that the principles I am using work 100%
correctly before making my code smarter.
Now, for a single processor the ghost atoms work the same as the periodic
images. That means if I wrap the molecules dimension I must get those of
ghost atoms (please correct me if I am wrong). I want to output the forces
for 300 atoms (including all atoms and their equivalent ghosts). I mean, if
I print (let's say) f[0][0], it's value show all the force counted for atom
*0* either considered as owned or ghost with other ids (just want to see
what the overall forces turned out to be)

i cannot answer your question. you are far too vague.
at the end of the force computation all required forces
are in the f- arrays of the "local" atoms. if newton_pair
is off, this is directly the case, if newton_pair
is on, then you need to wait until after
comm->reverse_comm() is called to have the contributions
computed on atoms in neighboring domains.

it is a *very* bad idea to think of lammps as a single processor code
that knows about the entire system.

axel.

5

Dear Axel,
Thanks for your guides. Turning the newton off I could match my forces with lammps. Now I am ready to implement the neighbor list.

However, doing so (turning the newton off) the pairwise energy changed from ~1063 to ~1120 kcal/mol. Is this a bug? Based on the documentation and your remarks everything should be the same with both newton on and off and only the method is less efficient, so what could be the reason?

Thanks again
Emily

6

Dear Axel,
Thanks for your guides. Turning the newton off I could match my forces with
lammps. Now I am ready to implement the neighbor list.

huh??? why do we need another neighbor list implementation?

However, doing so (turning the newton off) the pairwise energy changed from
~1063 to ~1120 kcal/mol. Is this a bug? Based on the documentation and your

changed where for what input?

remarks everything should be the same with both newton on and off and only
the method is less efficient, so what could be the reason?

the energy is the same unless there is a bug.
i am *very* tired of answering vague questions
about stuff that i don't ever get to see.

if you get the wrong results, most likely you
made a mistake. i am neither psychic nor
do i own a proper crystal ball.

if you want help, you have to help people
to help you by providing proper information.
if you don't want to do this, you have to
learn how to help yourself.

axel.