VMD-L Mailing List

From: Andrew Dalke (dalke_at_bioreason.com)
Date: Fri Mar 05 1999 - 14:54:07 CST

Axel Berg asked about:
> MSMS surface
> Does anyone know about this problem, and does anyone know a way to
> fix this?

  I don't have MSMS here so I cannot try it out, and it's been
several years since I last worked with it. I took a quick scan
of the code and it looks like the masking of atoms in done in
MSMS and not VMD (ie, for SURF, VMD writes the subset of atoms
to the SURF binary, which computes the whole surface of that
set).

  The code is in MSMSInterface.C, and there's some self-test
code still available in it. I see that it does a test where
some atoms are masked, and it checks that the right number of
atoms are detected:

| printf("Masking out some atoms\n");
| MSMSInterface msms4;
| int flgs4[] = {1, 0, 0, 1};
| if (msms4.compute(1.5, 1.0, n, NULL, xyzr, flgs4) > 0) {
| // need to see that everything in msms4 is in msms1

  but as you see, I use the first atom, so if there's a bug in MSMS
where the offset
> starts always from the first atom of the pdb file and not from
> the first atom of the particular selection.
then it wouldn't have been detected.

  You might try compiling MSMSInterface as a stand alone program
and try fiddling those flags to see what the results are. I can
offer no additional help.

> The second question I have concerns H-bonds. Does anyone know a
> way to write the H-bonds (partners), that VMD determines and shows
> on the screen, to a file?

  Alas, no. The hbond code is tied directly to the rendering code
(it's located in DrawMolItem.C). Also, the code isn't very rigourous;
it works pretty well for proteins, but in more general systems it
will put hbonds between obviously wrong places.

  That being said, the algorithm is pretty simple:

   for all selected non-hydrogens "i":
      for all other selected non-hydrogens "j" (j>i):
         if i is not bonded to j and distance(i, j) < maxdist:
            for all hydrogens "k" bonded to j:
               if angle (i,k,j) < maxangle:
                  draw hydrogen bond between i and j
            for all hydrogens "k" bonded to i:
               if angle (i,k,j) < maxangle:
                  draw hydrogen bond between i and j

and in theory the following VMD script should give the hbonds
between two selections (which in VMD C++ code is the same selection).
However, my testing shows this gives different results than
the C++ code, and I don't have the time to figure out exactly
what's the difference, so I'll leave it here as a starting point.

proc find_hbonds {sel1 sel2 maxdist maxangle} {
    if {[$sel1 molid] != [$sel2 molid]} {
        error "Must be from the same molecule"
    }
    # convert from degrees to radians
    global M_PI
    set maxangle [expr $maxangle / 180.0 * $M_PI]
    
    # get the bond list information for the two sets of atoms
    foreach atom [$sel1 get {index bondlist x y z}] {
        lassign $atom id blist x y z
        set bondlist($id) $blist
        set coord($id) "$x $y $z"
    }
    foreach atom [$sel2 get {index bondlist x y z}] {
        lassign $atom id blist x y z
        set bondlist($id) $blist
        set coord($id) "$x $y $z"
    }
    # find the hydrogens
    set hsel [atomselect [$sel1 molid] "hydrogen"]
    set hydrogens [$hsel list]
    # and get their coordinates
    foreach atom [$hsel get {index x y z}] {
        lassign $atom id x y z
        set coord($id) "$x $y $z"
    }
    
    # now do the search
    set atoms1 [$sel1 list]
    set atoms2 [$sel2 list]
    foreach i $atoms1 {
        foreach j $atoms2 {
            if {$i >= $j} {
                continue
            }
            # make sure they aren't bonded
            if {[lsearch $bondlist($i) $j] != -1} {
                continue
            }
            # check they are close enough
            if {[vecdist $coord($i) $coord($j)] >= $maxdist} {
                continue
            }
            # Find hydrogens bonded to j
            foreach other_atom $bondlist($j) {
                if {[lsearch $hydrogens $other_atom] == -1} {
                    continue
                }
                # it's a hydrogen, so find the angle
                set donortoH [vecnorm [vecsub $coord($other_atom) \
                                       $coord($j)]]
                set Htoacceptor [vecnorm [vecsub $coord($i) \
                                          $coord($other_atom)]]
                set ang [acos [vecdot $donortoH $Htoacceptor]]
                if {$ang > $maxangle} {
                    continue
                }
                puts "$j $other_atom $i"
                #draw line $coord($i) $coord($j)
            }
            # Find hydrogens bonded to i
            foreach other_atom $bondlist($i) {
                if {[lsearch $hydrogens $other_atom] == -1} {
                    continue
                }
                # it's a hydrogen, so find the angle
                set donortoH [vecnorm [vecsub $coord($other_atom) \
                                       $coord($i)]]
                set Htoacceptor [vecnorm [vecsub $coord($j) \
                                          $coord($other_atom)]]
                set ang [acos [vecdot $donortoH $Htoacceptor]]
                if {$ang > $maxangle} {
                    continue
                }
                puts "$i $other_atom $j"
                #draw line $coord($i) $coord($j)
            }
        }
    }
}

                                                Andrew