# special_bonds_dihedral_yes

Dear all: I'm using 'dihedral yes' to remove the effects of
special_bonds lj 0.0 0.0 0.0 on 1-4 neighbors that are not the
first-last atoms of any dihedral. However, in my case, the number of
1-4 neighbors becomes 0 after trimming the 1-4 neighbor list using
'dihedral yes'. This means 'special_bonds dihedral yes' does not find
any 1-4 neighbors that belong to a dihedral, thus "special_bonds lj
0.0 0.0 0.0" would affect all 1-4 neighbors.

The chain topology in my system is:

neighboring atoms along a chain are connected by one type of bonds,
i.e. linear topology,
atom i, i+1, i+2 are explicitly specified as atoms in an angle,
atom i, i+1, i+2, i+3 are explicitly specified as atoms in a dihedral,

then 1-5 neighbors (i, i+4) along the original chain are connected by
a different type of bonds.

Why does 'dihedral yes' fails to work with the original 1-5 neighbors

Best,
Ting

Dear all: I'm using 'dihedral yes' to remove the effects of
special_bonds lj 0.0 0.0 0.0 on 1-4 neighbors that are not the
first-last atoms of any dihedral. However, in my case, the number of
1-4 neighbors becomes 0 after trimming the 1-4 neighbor list using
'dihedral yes'. This means 'special_bonds dihedral yes' does not find
any 1-4 neighbors that belong to a dihedral, thus "special_bonds lj
0.0 0.0 0.0" would affect all 1-4 neighbors.

The chain topology in my system is:

neighboring atoms along a chain are connected by one type of bonds,
i.e. linear topology,
atom i, i+1, i+2 are explicitly specified as atoms in an angle,
atom i, i+1, i+2, i+3 are explicitly specified as atoms in a dihedral,

then 1-5 neighbors (i, i+4) along the original chain are connected by
a different type of bonds.

Why does 'dihedral yes' fails to work with the original 1-5 neighbors

i feel i sound like a broken record...

please provide a - preferably small - example
input that demonstrates the problem.

impossible to argue these things just hypothetically.
there are far too many things that go wrong.

axel.

Hi Axel: Thank you for your reply. The following is part of my input
script. It specifies the tables for different potential fields.

Dear T GE,

Dihedral keyword is affecting the 1-4 neighbors which are not defined in the system.

See the document below.

Thanks.

Special_bonds does not use the dihedral list to identify 1-4 neighbors
and reweight their interaction. Instead it finds 1-4 neighbors based
on the topology of bond connectivity. With dihedral yes, it removes
1-4 neighbors that are not first-last atoms of any dihedral from the
neighbor list. What remains should be a list of dihedral atoms, and
special_bonds would affect the interaction between them, as expected .

My concern is that in my case, if you look at the output file, the
number of 1-4 neighbors after the trim is zero. This means
special_bonds (special.cpp source file) does not find any 1-4 neighbor
that belongs to a dihedral. Therefore, I am worried that the
interaction within dihedrals is not reweighted. (In my case, the
normal lj interaction not deactivated.)

I would expect the number of 1-4 neighbors after the trim to be twice
the number of dihedrals. Just as the number of 1-3 neighbors after the
trim is twice the number of angles.

I tried to understand the protocol used to track 1-2, 1-3 and 1-4
neighbors in the special.cpp file. Not easy to follow...

Could someone please explain the way special.cpp trims the 1-4 neighbor list?

Best,
Ting

2013/1/22 Chang Woon Jang <[email protected]>:

T GE

I am sorry that I do not understand what you said.

With the LAMMPS document, I will explain what is going on your system.

From your data file, for example, you defined like

Dihedrals

1 3 1 2 3 4
2 1 2 3 4 5

For 1 to 5 atoms, if you did not define 2 1 2 3 4 5, pairwise interaction between atom 2 and atom 5 interaction should be affected with full weighting factor 1.0 because you used “dihedral yes”.

However, you did define this dihedral term 2 1 2 3 4 5. So, the pairwise interaction between atom 2 and atom 5 is affected by special_bonds lj 0.0 0.0 0.0.

As I understood from your data file, you defined all dihedrals from atom 1 to atom 4608 in sequence. Therefore, your “lj values (0.0 0.0 0.0)” made the neighbors 1-4 0 after dihedral trim (I guess the wording “after dihedral trim” means using and after “special_bonds” command). Your keyword “dihedral yes” did not affect your system.

I may be wrong. But with your data file and setting special_bonds in your input file, it makes sense to me based on the LAMMPS document.

He, Axel, will explain it anyway. I am also curious of my analysis of your output.

Thanks.

Changwoon Jang

Hi Axel: Thank you for your reply. The following is part of my input
script. It specifies the tables for different potential fields.

this is *useless* since it is incomplete.

axel.

T GE

I am sorry that I do not understand what you said.

agreed. it is very confusing and it doesn't make much sense.
but it seems to me the problem is really a badly designed
molecular model in combination with an incomplete
understanding of how exclusions are applied.

With the LAMMPS document, I will explain what is going on your system.

From your data file, for example, you defined like

Dihedrals

1 3 1 2 3 4
2 1 2 3 4 5
.....

For 1 to 5 atoms, if you did not define 2 1 2 3 4 5, pairwise interaction
between atom 2 and atom 5 interaction should be affected with full weighting
factor 1.0 because you used "dihedral yes".

However, you did define this dihedral term 2 1 2 3 4 5. So, the pairwise
interaction between atom 2 and atom 5 is affected by special_bonds lj 0.0
0.0 0.0.

this is correct. which means that both should be listed as "special dihedrals".

As I understood from your data file, you defined all dihedrals from atom 1
to atom 4608 in sequence. Therefore, your "lj values (0.0 0.0 0.0)" made the
neighbors 1-4 0 after dihedral trim (I guess the wording "after dihedral
trim" means using and after "special_bonds" command). Your keyword "dihedral
yes" did not affect your system.

this is not a correct analysis of the situation.
the real problem arises from the comment about
"adding bonds to model 1-5 neighbors".

from what i understood, there is not only a straight chain of bonds,
but *additional* bonds between atoms 1 and 5, 2 and 6, 3 and 7
and so on. *this* is the origin of the problem, because of that, the
bond topology is *massively* changed. the atoms 1,2,3,4,5 now also
form a ring and thus there are suddenly more topological angles
and dihedrals (e.g. 4,5,1 becomes a topological angle, and thus
the 1-4 interaction would be counted as both, a 1-3 and a 1-4 interaction.
the 1-3 takes precedence over the 1-4, so the dihedral 1,2,3,4 is
not even considered for a 1-4 interaction.

I may be wrong. But with your data file and setting special_bonds in your
input file, it makes sense to me based on the LAMMPS document.

yes, LAMMPS works correctly and as documented.
in order to demonstrate this, i am creating a set of
three data files with 3 different topologies made from
a linear chain of 11 atoms. a) has all topological dihedrals
defined, b) has one dihedral removed and c) has bonds
between 1-5 pairs added. the attached input script now
runs LAMMPS through the process of building the exclusions
for 3 special_bonds settings. default, dihedral yes and angle+dihedral yes

have a look and see.

axel.

p.s.: the attached files may also serve as an example,
how one can write minimal example inputs that demonstrate
how a feature does (not) work.
this was done my slightly modifying the VMD script code from here.

in.chain (2.23 KB)

data.chain-c (1.25 KB)

data.chain-b (1.19 KB)

data.chain-b (1.19 KB)

The following is a complete input script and the datafile, which
demonstrate that special_bonds doesn't work for the loop structure I
use. The structure is a single chain 1-2-3-4-5-6, with 1--5 and 2--6
connected. According to the screen output, 1-4 neighbors after the
trim is zero. My understanding is that this means special_bonds would
not affect the dihedrals I specified in the datafile.

The reason why I did this is that I need a pairwise potential between
1-5 neighbors that is different from the rest pairwise potential. What
I did before is modifing LAMMPS to calculate pairwise potentials based
on the atom IDs. If two neighbors belong to the same chain and are
1-5 neighbors, then they interact with a second pairwise potential.

I was wondering whether connecting 1-5 neighbors with tabulated bonds
(with large cutoff) can be another way to add special 1-5 pairwise
interaction. That's why I tried the wired loop structure. Seems that
the whole topology is complicated by this structure, and it doesn't
work.

Thanks,
Ting

T GE,

In your sample example, you defined pairwise cutoff distance 1.5 (e.g., lj/cut 1.5) but the pairwise distance between atom 1 and atom 4 is 3.

1 1 1 0.0 0.0 0.0
4 1 1 0.0 3.0 0.0

(I GUESS) That is why the interaction between atom 1 and atom 4 did not happen.

Thanks.

Changwoon Jang

The following is a complete input script and the datafile, which
demonstrate that special_bonds doesn't work for the loop structure I
use. The structure is a single chain 1-2-3-4-5-6, with 1--5 and 2--6
connected. According to the screen output, 1-4 neighbors after the
trim is zero. My understanding is that this means special_bonds would
not affect the dihedrals I specified in the datafile.

this is wrong. due to the additional bonds
you a) exclude the 1-5 non-bonded interaction
since it becomes a bond and excluded due
to a scaling factor of 0.0 and b) turn the 1-4
interaction of the dihedral 1-2-3-4 into a 1-3
interaction of the (topological) angle 1-5-4,
which is also excluded due to the scaling
factor of 0.0. the topological dihedrals remaining
before the trimming are *only* those including
the 1-5 bond (and 2-6) and since those are not
defined they all get removed from the list of
excluded interactions.

The reason why I did this is that I need a pairwise potential between
1-5 neighbors that is different from the rest pairwise potential. What
I did before is modifing LAMMPS to calculate pairwise potentials based
on the atom IDs. If two neighbors belong to the same chain and are
1-5 neighbors, then they interact with a second pairwise potential.

I was wondering whether connecting 1-5 neighbors with tabulated bonds
(with large cutoff) can be another way to add special 1-5 pairwise
interaction. That's why I tried the wired loop structure. Seems that
the whole topology is complicated by this structure, and it doesn't
work.

no, it *cannot* work. the problem is that you do this with a bond.
i can see multiple ways to support this, depending on how much
effort you want to invest and how clean a solution you want.

a) add a new interaction class for 1-5 interactions with support for
excluding them, too. that would be a *huge* amount of work
since you need to update the majority of files that are part of
lammps.

b) write a pair style that supports two parameter sets for a given
pair of atom types and you pick the one you want based on
both atoms having the same or a different molecule id.
this is a some additional work but probably the cleanest solution.

c) write a custom angle not bond style (so it doesn't add to the
bond topology) that then ignores the middle atom and implements
the modified pair interaction by applying the *difference* to the
regular non-bonded interaction for atoms of those types.
this is probably the least amount of programming, but
more work later.

d) find and modify a different code than lammps.
naaawwww....

axel.

Thank you, guys. Now I see why special_bonds works for the loop structure.

LAMMPS to support the feature in my system. With bonds connecting 1-5
neighbors, I just want to find another way of doing this.

You say "it *cannot* work. the problem is that you do this with a
bond." Is there anything that immediately excludes this option? I can
see the bond length could be a problem, but if I tabulate the bonding
potential and use a large cutoff, it still doesn't work?

Thanks,
Ting

Thank you, guys. Now I see why special_bonds works for the loop structure.

LAMMPS to support the feature in my system. With bonds connecting 1-5
neighbors, I just want to find another way of doing this.

again, you are not making sense here.

You say "it *cannot* work. the problem is that you do this with a
bond." Is there anything that immediately excludes this option? I can

i already explained in detail and with giving
examples why it cannot work... *twice*.
and gave you suggestions for workable
alternatives.

i am not going to explain it a third time.
i wouldn't know which language to use
without becoming offensive.

axel.

Okay. Thanks.