programming question - get bonds/angles/dihedrals for J-th atom

Dear All,

I have a question related to LAMMPS structure and programming,

I have a j-th atom, how could I find out in which bonds/angles/dihedrals it
belongs to?

in class Neighbor I see

The atom class stores this info for each atom.

num_bond[i] = # of bonds stored for (local) atom i
bond_type[i][j] = bond type for Jth bond of atom i
bond_atom[i][j] = global ID of Jth atom that atom i is bonded to

The caveat is that if newton bond is on (which it is by default),
then the I-J bond is only stored once, be either atom I of atom J.

If you set newton bond off, then if is stored twice, and angles
are stored 3 times, etc.

p/s/ is there any documentation (apart from notes in source code) for
LAMMPS,
like class tree, methods their description and general structure?

nope, sorry, but there is not

Steve

Dear Steve,

by "stored once" do you mean that
only one of 2 atoms will have a
non-zero num_bond[i] ?

Kind regards,
Denis.

exactly - this is so that the bond is only computed
once, and force imparted to both atoms I and J.

Steve

Also, there is an array int **special and int **nspecial (in atom.cpp)
that stores topology lists of all 1-2, 1-3, 1-4 neighbors. In these
lists J will be a 1st neighbor of I and I will also be a 1st neighbor
of J. So you could find "bonds" that way, if that is what you need.

Steve

Dear Steve,

Thank you for your answer.

I noticed that an atom ID = bond_atom[i][j] is global, as you said.
How do I convert it to a local, so that I can use that local index to access
xatom[ID_loc][k] and other fields? Probably there is a map from global to local.

I also went in debug mode through some code I have know,
and noticed that NeighList list_->numneigh could actually be zero
even for a simple example of two atoms, which are within lj/cut distance.
That field is accessed through a fix in AtC package after a 1 cycle Run command.

Is it possible that numneigh is zero even though two atoms are close enough?

I attach input and data files just for any case.

Thank you in advance,
Kind regards,
Denis.

in.2a (3.31 KB)

data.2a (365 Bytes)

Dear Steve,

To follow the question of neighbors:

Do I understand it correct that actually
no neighbors list is build for bonded atoms,
wherever they are coming from chains, like 1-2 + 2-3 + 3-4 or
from angle/dihedral potential like 1-2-3, 1-2-3-4.

it seems I have to stick to Atom's
int **nspecial;
int **special;
int maxspecial;

So, does special hold ID-s in the following way:

for 1-2: special[loc_ID][0 :n12 -1 ]
for 1-3: special[loc_ID][n12 :n123 -1 ]
for 1-4: special[loc_ID][n123 :n1234-1 ]

where

n12 = nspecial[loc_ID][0]
n123 = nspecial[loc_ID][1]
n1234 = nspecial[loc_ID][2]

?

If so, can I ask you please to briefly describe where should I get special weights.
I could only find an
int special_flag[4] (why 4?),
but not weights themselves.

Thank you in advance,

Kind regards,
Denis.

Comments below.

Steve

Dear Steve,

Thank you for your answer.

I noticed that an atom ID = bond_atom[i][j] is global, as you said.
How do I convert it to a local, so that I can use that local index to access
xatom[ID_loc][k] and other fields? Probably there is a map from global to
local.

Indeed - see atom->map().

I also went in debug mode through some code I have know,
and noticed that NeighList list_->numneigh could actually be zero
even for a simple example of two atoms, which are within lj/cut distance.
That field is accessed through a fix in AtC package after a 1 cycle Run
command.

The I-J pair is typically stored either with atom I or with atom J,
so that could be the reason.

Comments below.

Steve

Dear Steve,

To follow the question of neighbors:

Do I understand it correct that actually
no neighbors list is build for bonded atoms,
wherever they are coming from chains, like 1-2 + 2-3 + 3-4 or
from angle/dihedral potential like 1-2-3, 1-2-3-4.

it seems I have to stick to Atom's
int **nspecial;
int **special;
int maxspecial;

So, does special hold ID-s in the following way:

for 1-2: special[loc_ID][0 :n12 -1 ]
for 1-3: special[loc_ID][n12 :n123 -1 ]
for 1-4: special[loc_ID][n123 :n1234-1

where

n12 = nspecial[loc_ID][0]
n123 = nspecial[loc_ID][1]
n1234 = nspecial[loc_ID][2]

yes - that's the idea

?

If so, can I ask you please to briefly describe where should I get special
weights.
I could only find an
int special_flag[4] (why 4?),
but not weights themselves.

force.h has special_lj[4] and special_coul[4]. It is length 4,
b/c it is accessed out of the neighbor list with a 1,2,3 index.
Index 0 is not used.

Dear Steve,

Thank you for your answers.

During debugging i run into some things which are not clear to me:

My "toy" system consists of 5 atoms with 1-2 LJ bonds only: 1-2, 2-3, 3-4 (5th is not bonded).

I have two lists to build, one is full (from ATC fix) and another if half. They seem to be ordered correctly for list build (full first, then half).

Verlet::setup() calls for appropriate list build, the first one is Neighbor::full_bin()

As a result I have:
numneigh = {4,4,4,4,4}
firstneigh[0][0-2] = some arbitrary integers (apparently not initialized),

and only
firstneigh[0][3] = 4 -- which is exactly what i would expect : 1st atom has only 1 non-bonded neighbor - 5th atom (given that 5-th is inside LJ-cut+skin).

I don't know if full_bin() is expected to behave like that (contain some rubish numbers in firstneigh), but if yes, how do I get "real" non-bonded neighbors for a given local_ID ?

Kind regards,
Denis Davydov

Dear Steve,

I think I found out how to access atoms in a right way (looking at pair_lj_cut.cpp):

j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;

i assume these are conventions made for the whole code and
are actually true/correct for all interactions.

Kind regards,
Denis.

If you are using firstneigh and numneigh after a neighbor
list has been built, and you are using it like the pair routines
use it, then there is little chance that it is built incorrectly.
You must be interpreting its contents incorrectly.

Steve

The masking is for special_bonds flags that are stored
in hi-order bits (sometimes).

Steve

Yes, I agree.
I was wrong in interpreting its contents, assuming that
the list contains only non-bonded neighbors.

Kind regards,
Denis.