How to calculate the average bond angle

LAMMPS LAMMPS Mailing List Mirror
1

Hello,
I am trying to calculate the average bond angle in a metal-organic framework. I tried setting up a group (with atom types 3,5,7), hoping to calculate all of the bonds formed by atoms 5-3-7. However, atoms of type 5 are bonded to two atoms of type 3 (so I have a chain of atom types 3-5-3-7). Because I had specified the angle based on the group, LAMMPS was actually calculating all of the 3-5-3 angles and all of the 5-3-7 angles and averaging both sets together.

Here is the code I have been using for this:

group C1NZn type 3 5 7
…
compute C1NZnAngle C1NZn angle/local theta
fix C1NZnAnglePrint C1NZn ave/histo 1 100 100 100.0 140.0 20 c_C1NZnAngle mode vector file C1NZnAngles.txt ave running
compute C1NZnAngleAvg all reduce ave c_C1NZnAngle
fix 14 C1NZn ave/time {Nevery} {Nrepeat} ${Nfreq} c_C1NZnAngleAvg file C1NZnAngleAvg.txt

From the “fix C1NZnAnglePrint … “ command I was able to see that the bond distribution is bimodal, with peaks at the expected values of the 3-5-3 angles and the 5-3-7 angles.

Is there a way to calculate an average bond angle based only on the angle number? I do have the angles defined, so the 3-5-7 angles are angle type 6 (and the 3-5-3 angles are angle type 8).

Any insights into this problem would be greatly appreciated!

Marie

2

I don't think LAMMPS can make this distinction based on
angle type currently. I can
think of 2 ways it could be extended to allow this:

1) allow groups to be defined that only include atoms that
are part of angles of specified angle types. I'm thinking
this still might fail, if 3 atoms are all part of a 353 angle, but
all 3 are also part of other unwanted angles like a 335 ??

2) allow for an extra modifier on all the compute local
commands that work with angles (or dihedrals, etc)
that restricts the command to a subset of angle types

It might be better if there were some way to do sub-select
the values wanted when fix ave/histo uses them, but
I'm not seeing an easy way to do that.

Maybe other folks have suggestions - I need to think
about this some more.

Steve