### Author Topic: How can we get a bond angle distribution in ovito?  (Read 2641 times)

#### sbista

• Newbie
• Posts: 2
##### How can we get a bond angle distribution in ovito?
« on: February 16, 2018, 05:32:37 AM »
Hi,
I would like to know if we could write a python script to get a bond angle distribution for a selected particle types in ovito. I am interested in creating a histogram of bond angle distribution between certain particle types in my simulations.
Thanks

#### Alexander Stukowski

• Hero Member
• Posts: 638
##### Re: How can we get a bond angle distribution in ovito?
« Reply #1 on: February 16, 2018, 04:17:41 PM »
Hi there,

Here is a script I wrote some time ago for Ovito 2.9.0 to calculate the bond-angle distribution for a system (for demonstration purposes). It may serve you as a starting point to develop a solution for your specific problem. Disclaimer: I never used the script for real-word analyses and I never verified that it is working correctly. So use with caution.

-Alex

#### SC

• Newbie
• Posts: 7
##### Re: How can we get a bond angle distribution in ovito?
« Reply #2 on: October 24, 2018, 03:58:20 PM »
Dear Dr. Stukowski,

I used your script to calculate bond angle distribution. I get the following:

Bond angle cosine histogram:
[4423220       0       0       0      93    6950  147557  998011 2604010
2946143 1606165  457026   73560    6972     416      13       0       0
17     612   12790  140796  844062 2757451 4967918 5003619 2852902
945803  261040  284857  629331 1145734 1616595 1763957 1499768  995794
515809  211497   68083   17656    3582     580     560    4518   30945
162271  646447 1931808 4161468 6223576 6294421 4273684 1978276  647640
153969   27062    3385     308      14       3       0       0       0
0       0       0       0       0       0       8     689   22868
343535 2270218 6146446 6336389 2299107  264943    8670      72       0
0       0       0       0       0       5     979   55577 1012293
4272345 3194279  309097    2341       0       0       0       0       0
0]
Invalid Python script. It does not define the function modify().

However, I am not sure how this can be interpreted. Would you please give some explanation about this?

Bets wishes,

#### Constanze Kalcher

• Sr. Member
• Posts: 301
##### Re: How can we get a bond angle distribution in ovito?
« Reply #3 on: October 24, 2018, 04:50:13 PM »
Dear SC,

concerning the error message: did you just copy that code in a python script modifier in the GUI of OVITO? Note that this is a bash-script which is meant to be executed from the command line like this:
Code: [Select]
ovitos <name-of-script.py>
See the documentation about running scripts.

-Constanze
« Last Edit: October 24, 2018, 04:51:50 PM by Constanze Kalcher »

#### SC

• Newbie
• Posts: 7
##### Re: How can we get a bond angle distribution in ovito?
« Reply #4 on: February 17, 2019, 04:47:45 PM »
Dear Constanze,

Do you mean that I should run it in "ovitos" program? or Windows commend line?

I tried both but was not successful.

Best,

#### Alexander Stukowski

• Hero Member
• Posts: 638
##### Re: How can we get a bond angle distribution in ovito?
« Reply #5 on: February 17, 2019, 06:46:02 PM »
Hi SC,

Constanze went on vacation, so let me answer your question.

Yes, "ovitos" is a command line application for executing Python scripts that do run outside of the graphical Ovito user interface. It is typically used when automating analysis and rendering tasks which do not require user intervention.

On Windows, you can find "ovitos.exe" in the installation directory of Ovito. You should run it from the Windows command line prompt as specified by Constanze in order to see the console output your Python script produces.

-Alex

#### SC

• Newbie
• Posts: 7
##### Re: How can we get a bond angle distribution in ovito?
« Reply #6 on: February 17, 2019, 07:14:37 PM »
Dear Alex,

Thank you so much for your reply.

I have attached a file showing the results of the analysis. However, i am not sure how to interpret the results. I would like to get a graph like that of g(r) which can be obtained in GUI.

Best,

#### Alexander Stukowski

• Hero Member
• Posts: 638
##### Re: How can we get a bond angle distribution in ovito?
« Reply #7 on: February 19, 2019, 07:21:42 PM »
The script builds up a histogram of bond angles over the [-1,+1] value range of cos(theta). This histogram is stored in a Numpy array in the variable "angle_cosine_histogram". You can use the numpy.savetxt() function to write the numeric data to a text file in the current directory. For example, add the following line to the end of the script:

Code: [Select]
np.savetxt("distribution.txt", angle_cosine_histogram, delimiter=',')

You can then use your favorite plotting tool (e.g. Microsoft Excel, Gnuplot, etc.) to import this text file and create a graph from the numeric data. There is currently no direct way of exporting the visual plots you can find within the graphical version of Ovito as images.

#### stl4

• Newbie
• Posts: 4
##### Re: How can we get a bond angle distribution in ovito?
« Reply #8 on: March 11, 2019, 11:15:20 PM »
Hello,

I have a question regarding the script you had provided. If I am interested in calculating the distribution of bonding angles between Si-O-Si as well as O-Si-O.
1)Would I need to alter the following portion "create_bonds_mod.set_pairwise_cutoff('Si', 'O', 1.7)" to calculate the O-Si-O and Si-O-Si distribution? I assume that the previous line calculates both the O-Si-O and Si-O-Si angle since my histogram appears in the attached image.

2)If both angles are calculated, how can I separate those two angles since the lower end of Si-O-Si occurs at the upper end of O-Si-O?

Thank you for your assistance
« Last Edit: March 20, 2019, 09:55:22 AM by Constanze Kalcher »

#### stl4

• Newbie
• Posts: 4
##### Re: How can we get a bond angle distribution in ovito?
« Reply #9 on: March 19, 2019, 03:46:33 PM »
Hello Alexander,

I am still looking for some clarity regarding which angle is produced in the script that you have provided. Thank you for your time!

Sean

#### Constanze Kalcher

• Sr. Member
• Posts: 301
##### Re: How can we get a bond angle distribution in ovito?
« Reply #10 on: March 19, 2019, 06:43:13 PM »
Hello Sean,

the script loops over all particles regardless of their particle type and uses all bonds that you have created with the Create Bonds modifier.

1) yes, in this case you only need to create bonds between Si and O.

2) In order to obtain the two different bond angle distributions, you could try and modify the part of the script where you iterate over all particles
such that you keep track of two histograms, one which is incremented when the central particle is Si and one that is only incremented when
the central particle is O.

Code: [Select]
...
tprop = data.particles["Particle Type"]
Si_type_id = tprop.type_by_name('Si').id
O_type_id = tprop.type_by_name('O').id
for particle_index in range(data.number_of_particles):
...
if (tprop[particle_index] == Si_type_id):
angle_cosine_histogram_1 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]
elif (tprop[particle_index] == O_type_id):
angle_cosine_histogram_2 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]

-Constanze

#### vhponcev

• Newbie
• Posts: 1
##### Re: How can we get a bond angle distribution in ovito?
« Reply #11 on: April 18, 2019, 03:11:43 AM »
Dear Alexander,
I have been modifying your code for several, frames:
Code: [Select]
from ovito.io import *
from ovito.data import *
from ovito.modifiers import *
import numpy as np
import math
f1="../heating.dump"
node1 = import_file(f1, multiple_frames = True)
t=[]
#export frames separately
export_file(node1, "output.*", "lammps_dump", columns=
["Particle Identifier", "Particle Type", "Position.X", "Position.Y", "Position.Z"]
,multiple_frames=True)
# Number of bins of the histogram to create
bin_count = 1000
#array of angles
for j in range (bin_count):
k=(180/math.pi)*(math.acos(j/(bin_count/2)-1))
t.append(k)
angle_cosine_histogram_1 = np.zeros((bin_count,), int)
angle_cosine_histogram_2 = np.zeros((bin_count,), int)

#count in every exportd frame
for fram in range(node1.source.num_frames):
node = import_file("output."+str(fram), multiple_frames = True)

create_bonds_mod = CreateBondsModifier(mode=CreateBondsModifier.Mode.Pairwise)
create_bonds_mod.set_pairwise_cutoff('Type 2', 'Type 5', 2.3)
node.modifiers.append(create_bonds_mod)

# Compute normalized bond vectors
data = node.compute()
particle_positions = data.particle_properties.position.array
bonds_array = data.bonds.array
bond_vectors = particle_positions[bonds_array[:,1]] - particle_positions[bonds_array[:,0]]
bond_vectors += np.dot(data.cell.matrix[:,:3], data.bonds.pbc_vectors.T).T
bond_vectors /= np.linalg.norm(bond_vectors, axis=1).reshape(-1,1)

# Iterate over all particles and their bonds
bonds_enum = Bonds.Enumerator(data.bonds)
tprop = data.particle_properties["Particle Type"]
for particle_index in range(data.number_of_particles):

local_bonds = bond_vectors[np.fromiter(bonds_enum.bonds_of_particle(particle_index), np.intp)]
# Compute bond angle cosines of all bond pairs
angle_matrix = np.squeeze(np.inner(local_bonds[np.newaxis,:,:], local_bonds[:,np.newaxis,:]), axis=(0,3))

# Create linear array of bond angle cosines (all values in the upper half of matrix)
angle_cosines = angle_matrix[np.triu_indices(len(angle_matrix), k=1)]
#print (angle_matrix)

# Incrementally compute histogram of angle cosines for each center
if (tprop.array[particle_index] == 2):
angle_cosine_histogram_1 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]
elif (tprop.array[particle_index] == 5):
angle_cosine_histogram_2 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]
tot=angle_cosine_histogram_1+angle_cosine_histogram_2
result = np.column_stack((t,angle_cosine_histogram_1,angle_cosine_histogram_2,tot))
file2= open("out", 'w')
file2.write("theta(degree)  N(5-2-5)  N(2-5-2) "+'\n')
for j in range (bin_count):
file2.write(str(result[j][0])+" "+str(result[j][1])+" "+str(result[j][2])+" "+str(result[j][3])+" "+'\n')
file2.close()

The code above is working successfully, however, it is not efficiently, since I had to export each frame separately and then load those frames for calculating the angle distribution and added to the previous distribution, I would like to know another way to apply the code for several frames without exporting separately.

#### Constanze Kalcher

• Sr. Member
• Posts: 301
##### Re: How can we get a bond angle distribution in ovito?
« Reply #12 on: April 23, 2019, 10:44:40 AM »
Hi,

it does not make much sense to me why you import, export and then re-importing your files.
It's sufficient to just load your trajectory once and then loop over all the frames using the compute(<frame>) function which takes the current frame number as an argument.
Also you should only append the Create Bonds modifier to the pipeline once. Just do it before you loop over all frames, otherwise you add an additional Create Bonds modifier every time you go to a new frame.
Code: [Select]
node = import_file(f1,  multiple_frames = True)
#...
create_bonds_mod = CreateBondsModifier(mode=CreateBondsModifier.Mode.Pairwise)
create_bonds_mod.set_pairwise_cutoff('Type 2', 'Type 5', 2.3)
node.modifiers.append(create_bonds_mod)

for frame in range(node.source.num_frames):
data = node.compute(frame)
#....

-Constanze

#### francoMU

• Newbie
• Posts: 1
##### Re: How can we get a bond angle distribution in ovito?
« Reply #13 on: April 29, 2019, 01:26:24 PM »
In case somebody needs the script for the GUI and Ovito 3.x.x:

Code: [Select]
from ovito.data import *
import numpy as np
from numpy import linalg as LA
from matplotlib import pyplot as plt

def modify(frame, data):

plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.title('Bond Angle Distribution')

binwidth = 0.1
bins = np.arange(0,180,binwidth)
angle_hist = np.zeros(len(bins)-1,int)

topology = data.particles.bonds.topology
positions = data.particles.positions
bond_vectors = positions[topology[:,1]] - positions[topology[:,0]]

cell = data.cell
bond_vectors += np.dot(cell[:3,:3], data.particles.bonds.pbc_vectors.T).T

bond_vectors_norm = bond_vectors/np.expand_dims(np.linalg.norm(bond_vectors,axis=1),1)

bonds_enum = BondsEnumerator(data.particles.bonds)

for particle_index in range(data.particles.count):

index_bonds = np.fromiter(bonds_enum.bonds_of_particle(particle_index), int)

local_bonds = bond_vectors_norm[index_bonds]

angle_matrix = np.squeeze(np.inner(local_bonds[np.newaxis,:,:],
local_bonds[:,np.newaxis,:]), axis=(0,3))

angle_cosines = angle_matrix[np.triu_indices(len(angle_matrix), k=1)]

degree_angle = np.arccos(angle_cosines)*180/np.pi

angle_hist += np.histogram(degree_angle, bins=bins)[0]

c = np.column_stack((bins[:-1],angle_hist))
np.savetxt("file.txt",c)
plt.plot(bins[:-1],angle_hist)
plt.show()

#### Sacho

• Newbie
• Posts: 11
##### Re: How can we get a bond angle distribution in ovito?
« Reply #14 on: July 16, 2019, 09:31:13 AM »
Hello, can this program get a bond angle distribution centered on each atom? If I want to get the angle distribution of a specific bond pair1-2-3, firstly，create bonds between 1 and 2，between 2 and 3. Secondly，modify the part of the script where you iterate over all particles，code show as below：

#...
create_bonds_mod = CreateBondsModifier(mode=CreateBondsModifier.Mode.Pairwise)
create_bonds_mod.set_pairwise_cutoff('Type 1', 'Type 2', 2.3)
create_bonds_mod.set_pairwise_cutoff('Type 2', 'Type 3', 2.3)
node.modifiers.append(create_bonds_mod)

tprop = data.particles["Particle Type"]
for particle_index in range(data.number_of_particles):
...
if (tprop[particle_index] == 2_type_id):
angle_cosine_histogram_1 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]
else:
angle_cosine_histogram_2 += np.histogram(angle_cosines, bins=bin_count, range=(-1.0,1.0))[0]

Is this right?
Thank you very much.

#### Constanze Kalcher

• Sr. Member
• Posts: 301
##### Re: How can we get a bond angle distribution in ovito?
« Reply #15 on: July 16, 2019, 11:03:33 AM »
Hello,
Concerning your first question, yes, you can access individual bond vectors and calculate the angle between a pair of bond vectors. Here is the link to the scripting reference of the latest developer version of ovito:
https://ovito.org/manual/python/modules/ovito_data.html#ovito.data.Bonds

As for your second question, could you please again explain what you're trying to do?

-Constanze

#### Sacho

• Newbie
• Posts: 11
##### Re: How can we get a bond angle distribution in ovito?
« Reply #16 on: July 16, 2019, 11:26:49 AM »
Dear Constanze
I want to get the bond angle distribution centered on each atom, I saw your previous reply and want to ask you the code to solve this problem. The “firstly and secondly” in the last question is that I want to ask if the steps are correct. If it is right, I will use it.
Best wishes to you

#### Constanze Kalcher

• Sr. Member
• Posts: 301
##### Re: How can we get a bond angle distribution in ovito?
« Reply #17 on: July 16, 2019, 01:50:52 PM »
Sorry I have to ask again, what do you mean by "bond angle distribution centered on each atom"? Do you want to see how the bond angle between a specific bond pair is evolving with time (i.e you have more than just one frame)? Or do you want to have a histogram of the bond angle distribution created from all bond-vector pairs in your system (with the additional constraint that you want certain particle types as central particles)? In the latter case the code snippet you posted above should work.

-Constanze

#### Sacho

• Newbie
• Posts: 11
##### Re: How can we get a bond angle distribution in ovito?
« Reply #18 on: July 17, 2019, 03:25:44 AM »
Dear Constanze
Thanks very much. My promble has been solved, you are right, I want ro calculate the bong angle distribution for a specific bond pair. Your work is great.
-Sacho