# pizza.py

Dear Ben

I'm not sure it pizza.py can calculate bond-angles, (just because
you would have to indicate which atom triplets you want to use). But
if you can extract the relevant atom coordinates (using pizza.py or
other scripts), writing a program to calculate the bond-angles would
be an good programming exercise. (In case you don't have time for
that, read on.) There's an example how to do this for a linear chain
of atoms at the end of this email.

I wrote two relevant programs: coords2angles.py and dump2data.py

One of these programs is a tiny script which calculates bond-angles
between 3 atoms named "coords2angles.py". (There's also a
"coords2dihedrals.py".) You can download it here, and click on the
"dlpdb" link. ("coords2angles.py" is in the "extract_coords"
subdirectory)

http://moltemplate.org/other/index.html

----- documentation for coords2angles.py: ----

This program reads a 9-column numeric text file and prints a list
of bond-angles
The 9 numbers on each line represent the x,y,z coordinates of 3 atoms.
For each line this program computes the 3-body "bond-angle" between these atoms
and prints it to the standard out.
(When a line contains the wrong number of numbers, this program does not crash.
Instead an impossible value, -720, is printed.)

----- documentation for dump2data.py: ----

(Forgive me for forwarding this to the pizza.py user's list.)
Although pizza.py can extract atom coordinates, I'm too lazy to learn
pizza.py, so I'll tell you how to do everything using "dump2data.py"
which (like pizzapy), can extract atom coordinates, and save them as
either lammps DATA files or simple 3-column coordinate files ("RAW"
files). This program comes with moltemplate, so you should already
have it. If your LAMMPS dump file is named "traj.lammpstrj", you can
extract the coordinates this way:

dump2data.py -raw < traj.lammpstrj > traj.raw

This generates a 3-column text file containing the xyz coordinates on
each line of each atom (sorted by atomid).

data.py -raw < traj.lammpstrj > traj.raw

It assumes you are using "atom_style full". If you are using some
other atom style (eg "hybrid bond dipole"), then you can run it this
way:

dump2data.py -raw -atomstyle "hybrid bond dipole" \
< traj.lammpstrj > traj.raw

There is a blank line separating different frames in the trajectory,
so the number of lines in the file will be N*(n+1), where n is the
number of atoms, and N is the number of rames in the trajectory. If
this generates too many frames, you can limit this using:

data.py -raw -interval 10000 < traj.lammpstrj > traj.raw

to indicate the desired interval between frames (must be a multiple of
the save interval).

--- splitting text files and mergine lines ----

You will have multiple frames in your trajectory. You can use the
unix/linux "split" program:

http://linux.die.net/man/1/split

For example, if you have 1000 atoms, then (including the blank line
separator), you will have 1001 lines per frame. Split these using:
split -d -l 1001 - traj_frame < traj.raw

The program "coords2angles.py" expects 9 numbers per line, but the
coordinate files contain only 3 numbers per line (plus a blank line if
you use dump2data.py)

IF you want to calculate angles between successive atoms (in a linear
chain, for example), use this command to put the coordinates of 3
consecutive latoms onto one line.

awk '{if (NF==3) {l3=l2;l2=l1;l1=\$0} if (NR>2){print l3" "l2" "l1}}' \
< traj_frame00 > angles_frame00.dat

awk '{if (NF==3) {l3=l2;l2=l1;l1=\$0} if (NR>2){print l3" "l2" "l1}}' \
< traj_frame01 > angles_frame01.dat

:
etc

(Hopefully email did not split up those awk commands onto multiple
lines, or you get the idea anyway)

Hope this helps

Andrew

Andrew,
I have a dump file of x y and z coordinates and want to generate bond angle
distribution. Do you think this is possible with pizza.py?
Ben

Dear Ben

I'm not sure it pizza.py can calculate bond-angles

I feel like a moron.
I should google your question before boasting about my own programs.
(and probably so should you)
Yes pizza.py can do this. Check out:

http://pizza.sandia.gov/scripts.html

It looks like the only my software can do which pizza.py can not (yet)
do is calculate dihedral angles ("oords2dihedrals.py" in the
"extract_geometry" directory in the "dlpdb"), but that was not your
question (and it's trivial to modify the script which comes with
pizza.py to do this).

I hope I did not waste everybody's time.
Oops

Andrew