Welcome To Support Community

Materials Studio

Advanced Search
Ask Search:
Mohammad Abdul SalamMohammad Abdul Salam 

How to calculate total number of hydrogen bonds in a system

Hello Everyone,

I need to calculate the total number of hydrogen bonds in a system (periodic cell). My system contains different types of molecules. I have some system with two different molecules, some with three different molecules and some with just one type of molecules inside it. I don't need to know the hydrogen bonds between the two molecules, i need to know the total number of hydrogen bonds in these system to compare. Or else the total number of hydrogen bonds in the system over the period of time. I have trajectories run for 3000 ps with frame output at every 1000 steps and I am using Material Studio 7.0.

Hi Abdul
I haven't do that calculation, but I know that there is a script in (my version is 8.0)Accelrys\Materials Studio 8.0\share\examples\scripting\hbondstats.pl. I guess this maybe helpful

Mohammad Abdul SalamMohammad Abdul Salam
Hi Hanzi,

Thanks for the reply. I tried to run the script which you have suggested but it is getting failed. It is showing the error "Illegal division by zero at -e line 42". The -e line 42 is "$statsDoc->Cell($row, 1) = $totalLength/$row;". I am just copying and  pasting the script below, in this script I am just changing the -e line 17 (my $hbonds = $Documents{"urea.xsd"}->UnitCell->HydrogenBonds;) to my document name. 


# Purpose: Calculate statistics (Min/Max/Mean) on the Hydrogen Bonds
#            in a structure. Put the results in a Study Table.

use strict;
use warnings;
use MaterialsScript qw(:all);

#Initialize variables for the stats calculations
my $totalLength = 0;
my $minLength   = 99999.9; #Arbitrary Big Num
my $maxLength   = 0;
my $row = 0;

#Get all the HBonds in the UnitCell
my $hbonds = $Documents{"urea.xsd"}->UnitCell->HydrogenBonds;

#Create a new Study Table for the results
my $statsDoc = Documents->New("HBondStats.std");

foreach my $hbond (@$hbonds) {
    #Output the bond length for each HBond
    $statsDoc->Cell($row, 0) = "HBond $row";
    $statsDoc->Cell($row, 1) = $hbond->Length;

    #Update the statistics information
    $totalLength += $hbond->Length;

    if($hbond->Length < $minLength) {
        $minLength = $hbond->Length;
    elsif($hbond->Length > $maxLength) {
        $maxLength = $hbond->Length;


#printout the overall statistics
$statsDoc->Cell($row, 0) = "Average";
$statsDoc->Cell($row, 1) = $totalLength/$row;
$statsDoc->Cell($row, 2) = "Min";
$statsDoc->Cell($row, 3) = $minLength;
$statsDoc->Cell($row, 4) = "Max";
$statsDoc->Cell($row, 5) = $maxLength;

Jason DeJoannisJason DeJoannis
Hi Abdul,

I think that you don't have any hydrogen bonds, so it is attempting to divide by zero. Before running this script, you actually have to calculate the hydrogen bonds by clicking the Calculate Hydrogen Bonds button. The hydrogen bonds will be visible as blue dashed lines.User-added image
Jason DeJoannisJason DeJoannis
By the way, if you only want the number of hydrogen bonds, that is very easy:

printf "Nr hbonds %d\n", $doc->UnitCell->HydrogenBonds->Count;