angles cannot all be removed when bond is broken

Andrew:

Thanks very much for your kind help. I will try your suggestions.

The two ways you mentioned of changing the angle types, they both need to be done. One (anglelist[][]) is for angles in the neighbor list (temporal list), and the other (atom->angle_type[][]) is for the permanent angle list which will be written to the restart file.

If I understand right, GCMC involves insertions and deletions of molecules. When a new molecule is inserted, will there be new bonds and angles added to the permanent bond and angle list? If so, this actually resembles my problem. And one issues is, how the new bond and angle can be added on each processor if run parallely. My problem is, for an angle a1-a2-a3 (atom ids), if a1 and a2 are on processor P1 while a3 is on processor P2, when bond a1-a2 is broken the angle (a1-a2-a3) stored on P1 is set type 0 but the same angle stored on P2 is not be changed. If you have done some coding work regarding updating bond/angle information on multiple processors, I think that could shed some light on solving my current difficulty.

Best Regards
Shaorui

Andrew:

Thanks very much for your kind help. I will try your suggestions.

The two ways you mentioned of changing the angle types, they both need to be
done. One (anglelist[][]) is for angles in the neighbor list (temporal
list), and the other (atom->angle_type[][]) is for the permanent angle list
which will be written to the restart file.

If I understand right, GCMC involves insertions and deletions of molecules.
When a new molecule is inserted, will there be new bonds and angles added to
the permanent bond and angle list? If so, this actually resembles my
problem. And one issues is, how the new bond and angle can be added on each
processor if run parallely. My problem is, for an angle a1-a2-a3 (atom ids),
if a1 and a2 are on processor P1 while a3 is on processor P2, when bond
a1-a2 is broken the angle (a1-a2-a3) stored on P1 is set type 0 but the same
angle stored on P2 is not be changed.

After looking at your code, it occurs to me that you are not using any
MPI_ functions to communicate changes in the bonds to the other
processors. I would think that you should broadcast changes to the
bonds and angle-lists to the other processors using MPI_Bcast().

--- for example ---
For example, after you set an angle type to zero:
   atom->angle_type[ai1][am] = 0;
perhaps you should invoke MPI_Bcast() on it to share it with the other
processors.
   MPI_Bcast(&(atom->angle_type[ai1][am]),1,MPI_INT,0,world);
(I'm not sure if you need to broadcast any changes to "anglelist[][]"
to neighboring processors.)

Strangely (to me at least) I noticed that for some reason the current
version of BondQuartic::compute() in "bond_quartic.cpp" does not make
any MPI calls. Does the current version of "bond_style quartic" work
in parallel?

    I'm not sure if this will help, but if you have not done this yet,
take a look at the source code for "fix bond/break" and "fix
bond/create". Those files are at:
src/MC/fix_bond_break.cpp
src/MC/fix_bond_create.cpp
    (After a casual glance, I was not able to figure out what was going on.)

    The problem with your original post was that it's hard to ask
people to go through your code and look for errors. Perhaps once you
have a more tangible and easy-to-answer questions about something you
are confused about, you will get more replies.

   Alas, I suspect that more qualified people might not be replying to
your post because I did, which is a shame because I am not qualified
to talk about MPI.

Good luck!

Andrew

Strangely (to me at least) I noticed that for some reason the current
version of BondQuartic::compute() in "bond_quartic.cpp" does not make
any MPI calls. Does the current version of "bond_style quartic" work
in parallel?

yes, it works. and it doesn't need to use any communication,
since bond/quartic requires special 1 1 1

bond/handles excluded non-bonded interactions by
computing and subtracting them as part of the bond
(or not, in case the bond is broken).

the only information that would need to be communicated
would be the information for exclusions, but with 1 1 1,
there is no excluded interaction. anything else is handled
by the neighbor list code. if you have newton on, the bond
is only computed on one processors and if you have newton
off, the bond is broken on both processors. with a bond,
which involves only two atoms, this is all very simple.

add angle/dihedrals/impropers and exclusion factors != 1
to the mix and you have a big mess. now you'd *have*
to communicate any kind of information of topology changes
at the end of the bond compute. but how are you going
to change exclusions *after* you have applied them in pair?
you can do a similar trick as bond/quartic, but then you
have to also make more topology related changes and
then you have to keep track of what was changed when.
it hurts my brain even to ponder what complications
may arise with systems when the atoms are distributed
over more than two domains...

axel.

Thanks Axel. Although I did not understand every detail, I understood
enough to realize I am not ready to worry about messing with this
code, given my current knowledge of MPI and the magic of LAMMPS
neighboring code. Your email probably isn't very encouraging news to
Shaorui either.
It was interesting to learn about how bond style quartic works.

Andrew

Hi Andrew:

Thanks for your reply. Can I have a look your script (“gen_all_angles_topo.sh”)? I will keep working on modifying the quartic bond potential but in case I cannot get the job done, I will have to use some routine as you suggested…

Best Regards

Shaorui

Hi Axel:

In a previous posts by someone else out this issue, Steve suggested that modifying the angle type in the bond potential source file is probably the best option. But in case are a bond and its related angle are not stored on the same proc, for example, bond i1-i2 is stored on proc P0 and angle i1-i2-i3 is stored on proc P1, it seems impossible to modify the angle type if the bond i1-i2 is broken because, P1 can’t “see” i1-i2 and wouldn’t check if the angle should be modified. Am I right or I am just not aware of some way the job can be done?

Also, I was trying to learn something from the fix_bond_swap, since it involves change of angle partners after bond swap. Does it have communication among processors about updates of topology information?

Thanks.

Regards
Shaorui

Hi Axel:

In a previous posts by someone else out this issue, Steve suggested that
modifying the angle type in the bond potential source file is probably the
best option. But in case are a bond and its related angle are not stored on
the same proc, for example, bond i1-i2 is stored on proc P0 and angle
i1-i2-i3 is stored on proc P1, it seems impossible to modify the angle type
if the bond i1-i2 is broken because, P1 can't "see" i1-i2 and wouldn't check
if the angle should be modified. Am I right or I am just not aware of some
way the job can be done?

well, come to think of it. the best way to handle
the scenario you describe, might be to implement
this interaction entirely as an angle style and don't
even define the corresponding bonds (or set their
bond style to none) and then
- check if any of the bonds are broken
- if both, do nothing
- if one, compute the one bonded interaction and
  subtract the excluded 1-2 non-bonded force
- if both, handle both bonds as well as the angle
  potential and also subtract the 1-3 non-bonded force.

of course, this - same as bond/quartic - requires
special_bonds to be 1 1 1

but this way you have all bonds involved in an angle
on the same processor, and can also check which
components need to be computed/stored depending
on the newton settings.

Also, I was trying to learn something from the fix_bond_swap, since it
involves change of angle partners after bond swap. Does it have
communication among processors about updates of topology information?

no communication needed, since this is not called
during force computation but before reneighboring,
so the communication will come later.
please also note, that fix bond/swap requires that
1-3 non-bonded interactions are not excluded.

axel.

Hi.

I just sent you a version of the "gen_all_angles_topo.sh" script with
the most recent bug-fixes. (In a few days, I'll post it at
http://www.moltemplate.org/other.html)

(Changwoon was nice enough to test it. It worked for him, but it is
not fast. It is far too slow to run it every timestep. (Perhaps every
1000 timesteps?) It also requires that you define a consistent set of
rules which define ALL of the angles in your system. Not just the
angles involved in the quartic bonds.)

    To use it, you will have to also install "moltemplate" (moltemplate.org).

    Hopefully you can find another way.

Andrew

P.S.
(Until I update the web page, you will also need to copy the version
of "extract_lammps_data.py" from that email into moltemplate's "src"
directory. I'll update the web page later.)

    If you are desperate enough to consider using my script, then
perhaps you can somehow try getting your old customized version of
bond_quartic.cpp code to work by spamming redundant calls to ALL
processors after EVERY bond break using
   MPI_Bcast(&(atom->angle_type[ai1][am]),1,MPI_INT,0,world);
(or something similar). I suspect anything like that will still be
faster than periodically running my script, even if it's not as fast
as using the approach Axel and Steve suggested. (If it worked...)
Good luck.

Hello Axel:

I turned off the newton_bond (“newton on off” in the input script). It slowed down a little bit the simulation (about 10% more time for my system). The angles related to one bond can be all set to type 0 when that bond is broken. However, the simulation still got terminated with error “bond atom %d %d missing on proc %d at step” (or similar angle atom problem). The bond is one of those already broken at previous time steps. What could be the problem? Thanks.

Best
Shaorui

Hello Axel:

I turned off the newton_bond ("newton on off" in the input script). It
slowed down a little bit the simulation (about 10% more time for my system).
The angles related to one bond can be all set to type 0 when that bond is
broken. However, the simulation still got terminated with error "bond atom
%d %d missing on proc %d at step" (or similar angle atom problem). The bond
is one of those already broken at previous time steps. What could be the
problem? Thanks.

i don't know. i don't have a crystal ball or a sixth sense
and can't debug other peoples code without seeing them.

as i mentioned in my last e-mail, i believe you can avoid
the entire bookkeeping complications by implementing
your model as an angle style instead of a bond style.

turning off newton_bond won't solve any communication
problems. now you have to break the bond/angle consistently
on all involved processes.

axel.

Hello Axel:

In the LAMMPS source code, what is the criterion to invoke, this error (bond atoms
%d %d missing on proc %d at step %d)? Is it the bond length larger than pairwise
cutoff? Thanks.

Shaorui

Hello Axel:

In the LAMMPS source code, what is the criterion to invoke, this error (bond
atoms
%d %d missing on proc %d at step %d)? Is it the bond length larger than
pairwise
cutoff? Thanks.

no. the criterion is that atoms were not communicated at the re-neighboring.

axel.