Hi all: I am trying to modify LAMMPS codes to allow the use of cosine
angle potential with breakable quartic bond potential. Current version
of LAMMPS does not support broken bonds with 3 or 4 body interactions.

I notice that broken bonds are set as type 0 by LAMMPS, and that part
of bond potential would be skipped in later calculation. I am thinking
of using a similar method in angle potential calculation. If either
bond involved with an angle is of type 0, then this part of angle
potential is skipped. The question is can I refer to the bondlist
during angle potential calculation?

I have tried to expand the angle section in the input data file to
include the IDs of the two bonds associated with an angle. But then
realized that the bondlist might not be indexed in the original bond
ID. If the bondlist is indexed locally (I mean on individual
processors after ID maps), then using original bond IDs would make no
sense.

Another method I am thinking of using is altering atom types when
bonds are broken. This way since the anglelist already includes the
atom IDs, broken bonds can be identified in angle potential
calculation.

Does anyone have experience handling this before ? Any suggestions?
Thank you for your help.

I've recently implemented a fix "dynamic/angles". What it does is to determine periodicly which atoms form an angle and rebuilding the angle list accordingly. So for example if you have three types of angles in a system with three types of atoms (i.e. A (type 1),B (type 2) and C(type 3) and you want to have angle potentials for A-B-A triples (type 1), A-B-C triplets (type 2) and C-B-C triplets (type 3) you would define the following fix:

The first argument number is how often to rebuild the angle list, and the second number gives the number of angle types which will be managed by the fix (can be less than actual existing angle types) in this case three. For each angle type which should be managed you now give:
the type of the angle, the type of the central atom, the type of the first neighbor atom, the distance up to which the first neighbor is counted as neighbor, the type of the second neighbor atom and the distance up to which the second neighbor is counted as neighbor.

So the A-B-C angle definition is in the middle: 2 (type of angle) 2 (central atom type) 1 (first neighbor type) 1.5 (first neighbor max distance) 3 (second neighbor type) 1.2 (second neighbor max distance).

Maybe this is doing what you want?

Cheers
Christian

P.S. I have an additional options which allows to define different angle types depending on how the central atom is coordinated. I needed that for borate glasses where you trigonal BO3 units and tetragonal BO4 units.

The bond list is temporary (formed with neighbor lists).
The bond type stored by the atom is permanent and can
be accessed by the angle routines. Bond quartic changes
both, so you can access the permanent one.

A better way to do this however, is to have the bond routine
modify the angle type when it breaks a bond. Ditto for
diherdrals. That is non-trivial to do however.

Hi Steve: Thank you. I am thinking of using the following method. Have
the bond routine change the type of the two atoms in a broken bond
(permanently). Then in the angle potential calculation, if both atoms
in either bond are of the modified type, that angle is skipped. I will
preset the atom types and corresponding pairwise potential
coefficients based on my specific problem. Does this sound reasonable?
Anything needs special care?

You can do it however you wish. What you propose is
not a good general solution for LAMMPS since the atom
types are involved in many other parts of the code, e.g.
pair interactions. So what you are doing could have many
unintended side effects for other commands.

Hi Steve: Thank you. I see why changing the angle type is a better way
and actually necessary. Because skipping the part of calculation on an
angle potential without setting the angle type zero or negative would
lead to the error
"Angle atoms %d %d %d missing on proc %d at step".