Forum Navigation
You need to log in to create posts and topics.

Hydrogen Bonds

12

Within the Load Trajectory modifier settings you can either replace the 0 with a * character in the Search pattern field in the File sequence settings, or use the auto-generate button.

In case you have two different dump trajectories, i.e one with the particle position vectors and one with the dynamic bond information, the idea is then to use two Load trajectory modifiers in your pipeline to combine all this information.

Then in a next step, you can add a python script modifier to your pipeline. Click on edit script and open the script editor. You can take the script posted in this thread as a starting point for your own implementation of a bond counting function.
Let me know if you need further assistance with this.

Quote from Constanze Kalcher on August 21, 2020, 3:22 pm

Within the Load Trajectory modifier settings you can either replace the 0 with a * character in the Search pattern field in the File sequence settings, or use the auto-generate button.

In case you have two different dump trajectories, i.e one with the particle position vectors and one with the dynamic bond information, the idea is then to use two Load trajectory modifiers in your pipeline to combine all this information.

I tried it but it does not work in my Ovito GUI. It loads only one dump file "found 1 matching file". Maybe i need to update to ovito 3.2. I am using 3.0.

EDIT: I overcame this by saving the trajectories to a single dump file and loading it.

 

Quote from Constanze Kalcher on August 21, 2020, 3:55 pm

Then in a next step, you can add a python script modifier to your pipeline. Click on edit script and open the script editor. You can take the script posted in this thread as a starting point for your own implementation of a bond counting function.
Let me know if you need further assistance with this.

Can i write the whole procedure from the scratch by using only python scripting ? Start with reading the data file than add trajectory, create bonds and save the hydrogen bonds. It should be similar to the code you posted

from ovito.data import *
import numpy as np
from ovito.io import import_file, export_file
from ovito.modifiers import *

pipeline = import_file("/results/dump.*")

modifier = CreateBondsModifier()
modifier.set_pairwise_cutoff(3,11,3.5)
modifier.lower_cutoff = 1.1
pipeline.modifiers.append(modifier)

def hydrogen_bond(frame, data):

    #Get the particle type property
    particle_types = data.particles.particle_types
    #Get bonds topology and particle positions to compute all bond vectors present in your system
    topology = data.particles.bonds.topology
    positions = data.particles.positions
    bond_vectors = positions[topology[:,1]] - positions[topology[:,0]]
    bond_vectors += np.dot(data.cell[:3,:3], data.particles.bonds.pbc_vectors.T).T

    #Store this information as custom bond property, so they appear in the Data inspector
    data.particles_.bonds_.create_property("bond_vectors", data = bond_vectors)
    # Set up the bonds enumerator object.
    bonds_enum = BondsEnumerator(data.particles.bonds)
    # Counter
    bond_count = 0
    #for a,b in data.particles.bonds.topology:
    # print("Bond from particle %i to particle %i" % (a,b))
    
    # Loop over atoms.
    for particle_index in range(data.particles.count):
        # If the particle is a Hw atom continue and loop over bonds of current atom.
    
        if particle_types[particle_index] == 3:
            #Make sure that the current particle is connected
            bond_indices_list = [bond_index for bond_index in bonds_enum.bonds_of_particle(particle_index)]

            if len(bond_indices_list) >= 1:

                neighbor_particle_indices = []
                for bond_index in bonds_enum.bonds_of_particle(particle_index):

                # Obtain the indices of the two particles connected by the bond:
                a = topology[bond_index, 0]
                b = topology[bond_index, 1]

                if a == particle_index:
                    neighbor_particle_indices.append(b)

                else:
                    neighbor_particle_indices.append(a)

                #Verify that the particle types of the 2 neighboring atoms are actually to type Ow

                if( (particle_types[neighbor_particle_indices] == np.array([11,11])).all() ):
                    #Calculate angle
                    v_1 = bond_vectors[bond_indices_list[0]]
                    v_2 = bond_vectors[bond_indices_list[1]]
                    angle = np.arccos(np.clip(np.dot(v_1, v_2)/np.linalg.norm(v_1)/np.linalg.norm(v_2), -1.0,1.0 ))
                    
                    if np.degrees(angle) >= 120 :
                        bond_count+=1
    print("# Ow-Hw-Ow bonds:{}".format(bond_count))
    #Store result as global attribute so it appears in the Data inspector in the GUI
    data.attributes["Ow-Hw-Ow"] = bond_count

pipeline.modifiers.append(hydrogen_bond)

export_file(pipeline, "results/hb.txt",
format = "txt/attr",
columns = ["Timestep", "Ow-Hw-Ow"],
multiple_frames = True)

 

I tried like this but i am getting bonds not only between Hw and OW (3 and 11). How can i set CreateBondsModifier.Mode.Pairwise instead of Uniform ?

Hi,

ok great, I can see you managed to set up a stand-alone batch script. Note that you can change the mode of operation like this:

modifier = CreateBondsModifier()
modifier.mode = CreateBondsModifier.Mode.Pairwise

https://www.ovito.org/docs/current/python/modules/ovito_modifiers.php?highlight=create%20bonds#ovito.modifiers.CreateBondsModifier

Quote from Constanze Kalcher on August 24, 2020, 1:39 pm

Hi,

ok great, I can see you managed to set up a stand-alone batch script. Note that you can change the mode of operation like this:

modifier = CreateBondsModifier()
modifier.mode = CreateBondsModifier.Mode.Pairwise

https://www.ovito.org/docs/current/python/modules/ovito_modifiers.php?highlight=create%20bonds#ovito.modifiers.CreateBondsModifier

Thank you. Now i am getting my bonds between hydrogen and oxygen (3 and 11) and if i have multiple bonds how can i calculate the angle if i have e.g. 3 bonds or just one bond ?

PS: i do not have the insert code function next to insert/edit link. It is the insert/edit image function.

 

The Suppress inter-molecular bonds option means i have no bonds between different molecules but for water i need to exclude bonds within the molecule and include bonds between different water molecules.

if i have multiple bonds how can i calculate the angle if i have e.g. 3 bonds or just one bond ?

Correct me if I'm misunderstanding your question, but you can of course only calculate an angle between 2 bond vectors at a time, so you should make sure that the bond list of each particle in question has more than one entry. This is already included in the script you're using.

If the bond list of an atom contains more than 2 entries and you would like to calculate all angles between all possible combinations of bond pairs, you could extend your script after line 57.

PS: i do not have the insert code function next to insert/edit link. It is the insert/edit image function.

Thanks for letting me know about this, usually your post editor window should give you the code formatting options marked in the attached screenshot.

The Suppress inter-molecular bonds option means i have no bonds between different molecules but for water i need to exclude bonds within the molecule and include bonds between different water molecules.

If you have worked with molecule identifiers, then you could use those to verify if a bond connects two particles that are not part of the same molecule. Since you know the indices of the particles that are connected by each bond, you can use this information to look up the corresponding Molecule Identifier property entries. Otherwise, your only option is to use the lower-cutoff option in the Create bonds modifier as discussed above.

Uploaded files:
  • code-formatting.png

from ovito.data import *
import numpy as np
from ovito.io import import_file, export_file
from ovito.modifiers import *

pipeline = import_file()
traj_mod = LoadTrajectoryModifier()
traj_mod.source.load()
pipeline.modifiers.append(traj_mod)

modifier = CreateBondsModifier()
modifier.mode = CreateBondsModifier.Mode.Pairwise
modifier.set_pairwise_cutoff(3,11,3.5)
modifier.lower_cutoff = 0.0
pipeline.modifiers.append(modifier)

def unit_vector(vector):
""" Returns the unit vector of the vector. """
return vector / np.linalg.norm(vector)

def angle_between(v1, v2):
""" Returns the angle in radians between vectors 'v1' and 'v2'::

>>> angle_between((1, 0, 0), (0, 1, 0))
1.5707963267948966
>>> angle_between((1, 0, 0), (1, 0, 0))
0.0
>>> angle_between((1, 0, 0), (-1, 0, 0))
3.141592653589793
"""
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

#file = open("neighbo_list.txt", "w")
def hydrogen_bond(frame, data):

#Get the particle type property
particle_types = data.particles.particle_types
#Get the particle molecule id
molecule_id = data.particles['Molecule Identifier']
#Get bonds topology and particle positions to compute all bond vectors present in your system
topology = data.particles.bonds.topology
positions = data.particles.positions
bond_vectors = positions[topology[:,1]] - positions[topology[:,0]]
bond_vectors += np.dot(data.cell[:3,:3], data.particles.bonds.pbc_vectors.T).T
#Store this information as custom bond property, so they appear in the Data inspector
data.particles_.bonds_.create_property("bond_vectors", data = bond_vectors)
# Set up the bonds enumerator object.
bonds_enum = BondsEnumerator(data.particles.bonds)
# Counter
bond_count = 0
#for a,b in data.particles.bonds.topology:
# print("Bond from particle %i to particle %i" % (a,b))

# Loop over atoms.
for particle_index in range(data.particles.count):
# If the particle is a Hw atom continue and loop over bonds of current atom.
if particle_types[particle_index] == 3:
#Make sure that the current particle is connected
bond_indices_list = [bond_index for bond_index in bonds_enum.bonds_of_particle(particle_index)]
if len(bond_indices_list) >= 1:

neighbor_particle_indices = []
for bond_index in bonds_enum.bonds_of_particle(particle_index):
# Obtain the indices of the two particles connected by the bond:
a = topology[bond_index, 0]
b = topology[bond_index, 1]
if a == particle_index:
neighbor_particle_indices.append(b)
else:
neighbor_particle_indices.append(a)

#Get covalent bond to Ow of same water Molecule (Ow-Hw)

for i in range(len(neighbor_particle_indices)):
if molecule_id[neighbor_particle_indices[i]] == molecule_id[particle_index]:
cov_bond_index = bond_indices_list[i]
cov_oxygen = neighbor_particle_indices[i]
cov_bond_vector = bond_vectors[bond_indices_list[i]]

#Delete covalent bond from bond_indices_list
bond_indices_list.remove(cov_bond_index)
#Delete covalent Ow from neighbor_particle_indices
neighbor_particle_indices.remove(cov_oxygen)
#Verify that the particle types of the 2 neighboring atoms are actually to type Ow

for i in range(len(bond_indices_list)):
#Calculate angle between Donator-Hw and Hw-Acceptor
hb_vector = bond_vectors[bond_indices_list[i]]
angle = angle_between(cov_bond_vector, hb_vector)

if np.degrees(angle) >= 120:
bond_count+=1

# if( (particle_types[neighbor_particle_indices] == np.array([11,11])).all() ):
# #Calculate angle
# v_1 = bond_vectors[bond_indices_list[0]]
# v_2 = bond_vectors[bond_indices_list[1]]
# angle = np.arccos(np.clip(np.dot(v_1, v_2)/np.linalg.norm(v_1)/np.linalg.norm(v_2), -1.0,1.0 ))

# if np.degrees(angle) >= 120 :
# bond_count+=1
#print("# Ow-Hw-Ow bonds:{}".format(bond_count))
#Store result as global attribute so it appears in the Data inspector
data.attributes["Ow-Hw-Ow"] = bond_count
# data.attributes["Hb"] = neighbor_particle_indices

pipeline.modifiers.append(hydrogen_bond)

export_file(pipeline, "",
format = "txt/attr",
columns = ["Timestep", "Ow-Hw-Ow",],
multiple_frames = True)

Seems like this is working for me. Is there a possibility to write out for each timestep the number of hydrogen bonds for each oxygen atom ?

PS: That Option for inserting the code is not available for me

I would like to draw your attention to the new Bond Analysis modifier,

https://www.ovito.org/docs/current/particles.modifiers.bond_analysis.php and especially the option "Partitioned distributions -> Particle Type" which will allow you to easily obtain bond angle histograms for each combination of three particle types and bond length histograms for each pair-wise combination of particle types.

12