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

Hydrogen Bonds

12

Hello,

I am wondering if there is a way to calculate number of Hydrogen bonds (and their details) using Ovito, imposing two conditions, the distance between donor and receiving atoms is less than 3A and the angle is more than 120 degrees.

I have a hydroxylated Silica System with water. So there would be three types of hydrogen bonds, one within water (Ow-Hw-Ow), Silica walls with water (O-H-Ow) and Water with Silica walls (Ow-Hw-O).

(Silica Atom Type: Si-1, O-2, H-3, Water Atom Type: Ow-4, Hw-5)

Is there a way to calculate all three types individually.

Any help/suggestions are welcome and greatly appreciated.

Thank you

Pranay

Hello Pranay,

this is possible but it would require quite some python coding on your part. The idea would be to first create Bonds using the CreateBonds Modifier, see https://www.ovito.org/docs/current/particles.modifiers.create_bonds.php. Simply activate the "Pair wise cutoffs" option and edit only those cutoff values in the table for pairs you would like to see bonded.

Then in a next step, add a Python script modifier to your pipeline, https://www.ovito.org/docs/current/particles.modifiers.python_script.php. Please also have a look at the introduction on how to write user-defined modifier functions  https://www.ovito.org/docs/current/python/introduction/custom_modifiers.php for more details on the topic.

Let's start with the Ow-Hw-Ow bond-type: What you can do now, is to loop over all particles of type Hw and look up their bonds. This can be conveniently done using a BondsEnumerator and it's bonds_of_particle() method.
https://www.ovito.org/docs/current/python/modules/ovito_data.php#ovito.data.BondsEnumerator

If you find two bonds, you need to verify that the atoms connected by these bonds are actually of type Ow and Hw. If so you continue, and since you also said you only want to take into account bonds with a bond angle larger than 120 degrees, you then need to calculate the angle between those two bonds vectors. How to compute the bond vectors is explained here: https://www.ovito.org/docs/current/python/modules/ovito_data.php#ovito.data.BondsEnumerator

Finally, if the computed angle is larger than 120 degrees, you should count that bond pair.

Here is an example of such a modifier, if you like you can take this example, try it in the GUI and continue developing it.

from ovito.data import *
import numpy as np

def modify(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
 
    # 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] == 5:
            #Make sure that the current particle is connected to 2 bonds 
            bond_indices_list = [bond_index for bond_index in bonds_enum.bonds_of_particle(particle_index)]
            if len(bond_indices_list) == 2:
                                    
                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([4,4])).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

Hope that helps and let me know if you have questions or suggestions,

-Constanze

By the way, if you only want to keep those bonds in the GUI that fulfill your criteria, you should keep track of the bond_indices of all "valid bonds" and then later create a selection based on this information and delete all other bonds.

Hi Constanze,

Thank you very much for providing such useful information. I have two questions:

  1. Based on the code, does it mean that we have to manually calculate the maximum pair-wise cutoff for Ow and Hw according to the given distance and angle constraints in order to create and show all of the possible 'H Bond' in Ovito? Specifically, the Ow-Hw bond length within the water is about 1 angstrom. Considering the constraints of 3-angstrom distance and 120-degree angle, we should set about 2.4 angstroms for the pair-wise cutoff of Ow and Hw. I am not sure if I understand this procedure correctly.
  2. Code Line 27: 'Make sure that the current particle is connected to 2 bonds': Actually, I think a Hw atom of a water molecule could be bonded to more than one Ow atom of other water molecules (after create bond modifier). It seems that in such a case, the code excludes this Hw atom and the corresponding possible formed H bond involving this Hw atom because it is connected to more than 2 bonds. Should we take this kind of situation into account?

Many thanks in advance for any reply or comments.

 

Regards,

Lance

Hi Lance,

The way that I understood Pranay’s question, was that he wanted to know how to count the number of bonds (for different bond pair combinations), that exhibit a bond length < 3 A  and form an angle > 120 degrees and the answer above gives an example an how to do that for all water molecules.

I completely agree with you, that great care should be taken when generating these bonds with OVITO (if this information is not already part of your input data). You as a user need to know at what cutoff distance you consider two atoms bonded or not and should set the pair-wise cutoffs in the Create Bonds Modifier according to the problem you are studying.

In order to calculate an angle between to bond vectors, you should make sure that the center particle is connected to more than 1 other atom. If you don't want to restrain your analysis to particles that are bonded to two particles only, you could change line 29 to

if len(bond_indices_list) > 1:

and add another loop to iterate over all possible bond pair combinations.

Hope that helps,

Constanze

Quote from Constanze Kalcher on August 10, 2020, 12:52 pm

Hi Lance,

The way that I understood Pranay’s question, was that he wanted to know how to count the number of bonds (for different bond pair combinations), that exhibit a bond length < 3 A  and form an angle > 120 degrees and the answer above gives an example an how to do that for all water molecules.

I completely agree with you, that great care should be taken when generating these bonds with OVITO (if this information is not already part of your input data). You as a user need to know at what cutoff distance you consider two atoms bonded or not and should set the pair-wise cutoffs in the Create Bonds Modifier according to the problem you are studying.

In order to calculate an angle between to bond vectors, you should make sure that the center particle is connected to more than 1 other atom. If you don't want to restrain your analysis to particles that are bonded to two particles only, you could change line 29 to

if len(bond_indices_list) > 1:

and add another loop to iterate over all possible bond pair combinations.

Hope that helps,

Constanze

How can i exclude the already bonded water molecules (Hw to Ow) ? Is this even possible or just by setting the lower cutoff to e.g. 1.1 ?

Hi Betim,

could you give me some more explanation about what you would like to do?

Quote from Constanze Kalcher on August 21, 2020, 12:27 pm

Hi Betim,

could you give me some more explanation about what you would like to do?

i am trying to do the whole hydrogen bonding procedure with python scripting. If i do create bonds now within 3.5 angstrom i will have Hw-Ow bonds within the water molecules (Hw-Ow-Hw)  which are covalent bonds and not hydrogen bonds. Is it possible to exclude this and include only hydrogen bonds of hydrogens with oxygens of other molecules ?

Okay thank you, I see. Your initial suggestion of using the lower cutoff is probably a good idea. But without knowing more about your input data, it's difficult to give you any more specific instructions. In case you have used molecule identifiers, you could make use of the option "Suppress inter-molecular bonds" in the Create bonds modifier, e.g.

One more question, though. Could you tell me a little bit more about your input data? Have you used lammps for your simulation? If so, you might find the section "Varying bond connectivity" in the documentation of the Load Trajectory modifier helpful.

Quote from Constanze Kalcher on August 21, 2020, 2:44 pm

Okay thank you, I see. Your initial suggestion of using the lower cutoff is probably a good idea. But without knowing more about your input data, it's difficult to give you any more specific instructions.

One more question, though. Could you tell me a little bit more about your input data? Have you used lammps for your simulation? If so, you might find the section "Varying bond connectivity" in the documentation of the Load Trajectory modifier helpful.

I am using lammps and now i loaded my data file and how can i include now all my dump files ? The gui is including only my first "dump.0" file with the load trajectory modifier ?

Now i created bonds and did the "supress intermolecular bonding" for my Hw-Ow hydrogen bonding. How is it possible now to translate this into a python script and write the results out ? Is it similar to the script you put in ? Thank you for the help.

12