Scripts/formula.rb

   1 #
   2 #  main.rb
   3 #
   4 #  Created by Toshi Nagata.
   5 #  Copyright 2008 Toshi Nagata. All rights reserved.
   6 #
   7 # This program is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation version 2 of the License.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 class Molecule
  17
  18   @@known_fragments = Hash.new
  19   @@known_fragment_names = []
  20
  21   #  The instance variable @dummies is used to retain the indices of the detachable
  22   #  atoms in the fragment. For example, a methyl group "CH3" is internally
  23   #  represented by a 5-atom molecule "Du-CH3" where Du is a dummy atom. When
  24   #  a molecule is returned from the method Molecule.from_formula, the dummy atom
  25   #  is replaced with a hydrogen atom but the position of the dummy atom is retained
  26   #  as @dummies[0].
  27
  28   attr_accessor :dummies
  29
  30   #  Find dummy atoms; the atoms whose element is "Du" and name begines with "_"
  31   #  are regarded as dummy atoms.
  32   def find_dummy_atoms(how_many = nil)
  33     d = []
  34     atoms.each { |ap|
  35           if ap.element == "Du" && ap.name[0] == ?_
  36             d.push(ap.index)
  37                 if how_many != nil && d.length >= how_many
  38                   break
  39                 end
  40           end
  41         }
  42         d
  43   end
  44
  45   #  Returns the first atom connected to the specified atom.
  46   def root_atom(dummy_atom)
  47     return atoms[dummy_atom].connects[0]
  48   end
  49
  50   #  Register a fragment in @@known_fragment and rebuild @@known_fragment_names
  51   def Molecule.register_fragment(name, fragment)
  52     @@known_fragments[name] = fragment
  53         idx = @@known_fragment_names.length
  54         @@known_fragment_names.each_with_index { |n,i|
  55           if n.length < name.length
  56             idx = i
  57                 break
  58           end
  59         }
  60         @@known_fragment_names[idx, 0] = name
  61   end
  62
  63   def Molecule.known_fragment(name)
  64         @@known_fragments[name].dup
  65   end
  66
  67   def Molecule.lookup_fragment_from_string(str)
  68         @@known_fragment_names.each { |name|
  69           if str.rindex(name, 0) == 0
  70                 return name
  71           end
  72         }
  73         nil
  74   end
  75
  76   #  Create a fragment containing an atom and appropriate number of dummy atoms.
  77   #  The argument valence specifies the number of dummy atoms. Length is the
  78   #  bond length, and degree is the angle Du-X-Du (given in degree, not radian).
  79   def Molecule.atom_fragment(name, valence, length = nil, degree = nil)
  80     fragment = Molecule.new
  81         fragment.add_atom(name, "", name)  #  Atom type will be guessed later
  82         length = 1.0 if length == nil
  83         if valence > 0
  84           fragment.add_atom("_1", "", "Du", length, 0)
  85           if valence > 1
  86             if degree == nil
  87                   degree = 180.0 if valence == 2
  88                   degree = 120.0 if valence == 3
  89                   degree = 109.2 if valence == 4
  90                 end
  91                 fragment.add_atom("_2", "", "Du", length, 0, degree, 1)
  92                 if valence > 2
  93                   if valence == 3
  94                     sn = Math.sin(degree * Deg2Rad)
  95                         cs = Math.cos(degree * Deg2Rad)
  96                         acs = cs * (1 - cs) / (sn * sn)
  97                         if (acs < -1 || acs > 1)
  98                           dihed = 180.0
  99                         else
 100                           dihed = Math.acos(acs) * Rad2Deg
 101                         end
 102                   elsif valence == 4
 103                     dihed = 120.0
 104                   end
 105                   fragment.add_atom("_3", "", "Du", length, 0, degree, 1, dihed, 2)
 106                   if valence > 3
 107                     fragment.add_atom("_4", "", "Du", length, 0, degree, 1, -dihed, 2)
 108                   end
 109                 end
 110           end
 111         end
 112         fragment.dummies = fragment.find_dummy_atoms
 113         fragment
 114   end
 115
 116   #  Examine whether the given atom is a carbonyl carbon
 117   def is_carbonyl(idx)
 118     if atoms[idx].element == "C"
 119           atoms[idx].connects.each { |idx2|
 120             if atoms[idx2].element == "O" && atoms[idx2].connects.length == 1
 121                   return true
 122                 end
 123           }
 124         end
 125         return false
 126   end
 127
 128   def guess_bond_length(e1, e2, mult)
 129     if e1 > e2
 130           ee = e1
 131           e1 = e2
 132           e2 = ee
 133         end
 134         if (e1 == "C" && e2 == "C")
 135           len = (mult == 3 ? 1.18 : (mult == 2 ? 1.39 : 1.54))
 136         elsif (e1 == "C" && e2 == "H")
 137           len = 1.08
 138         elsif (e1 == "C" && e2 == "O")
 139           len = (mult == 2 ? 1.23 : 1.41)
 140         elsif (e1 == "C" && e2 == "N")
 141           len = (mult == 2 ? 1.29 : 1.48)
 142         elsif (e1 == "H" && e2 == "O")
 143           len = 0.96
 144         elsif (e1 == "H" && e2 == "N")
 145           len = 1.01
 146         elsif (e1 == "C" && e2 == "Cl")
 147           len = 1.79
 148         elsif (e1 == "C" && e2 == "S")
 149           len = 1.82
 150         elsif (e1 == "O" && e2 == "S")
 151           len = (mult == 2 ? 1.60 : 1.80)
 152         else
 153           len = 1.5
 154         end
 155         return len
 156   end
 157
 158   #  Make the 'root' atom trigonal (sp2). The atoms in the dummies array are
 159   #  considered as one dummy atom.
 160   def make_sp2(root, dummies, vec)
 161     nd = atoms[root].connects.select { |n| !dummies.include?(n) }
 162         return if nd.length != 2
 163         r0 = atoms[root].r
 164         r1 = atoms[nd[0]].r - r0
 165         ax = r1.cross(atoms[nd[1]].r - r0)
 166         if (ax.length < 1e-4)
 167           ax = r1.cross(atoms[dummies[0]].r - r0)
 168           if (ax.length < 1e-4)
 169             ax = r1.cross(Vector3D[1, 0, 0])
 170                 if (ax.length < 1e-4)
 171                   ax = r1.cross(Vector3D[0, 1, 0])
 172                   if (ax.length < 1e-4)
 173                     raise("Cannot determine axis for rotating fragment")
 174                   end
 175                 end
 176           end
 177         end
 178         r2 = atoms[nd[1]].r - r0
 179         cs = r1.dot(r2) / (r1.length * r2.length)
 180         ang = Math.acos(cs) * Rad2Deg
 181         rotate(ax, 120.0 - ang, r0, fragment(nd[0], root))
 182         v = ((atoms[nd[0]].r - r0).normalize + (atoms[nd[1]].r - r0).normalize) * (-1.0)
 183         vec.x = v.x
 184         vec.y = v.y
 185         vec.z = v.z
 186         return vec
 187   end
 188
 189   #  Add a molecular fragment to self. The new connection is created at the position
 190   #  of the dummy atoms. The 'first' dummy atom in self is kept untouched wherever possible.
 191   #  Mult (multiplicity) is the number of dummy atoms removed from each fragment.
 192   #  If unspecified, it is the same as the number of the dummy atoms in the added fragment.
 193   def add_formula(mol, mult = nil)
 194     natoms_s = self.natoms
 195         natoms_m = mol.natoms
 196         #  Find list of dummy atoms
 197     md = mol.find_dummy_atoms
 198         len = md.length
 199         if len == 0
 200           raise "The given molecule has no dummy atoms"
 201         end
 202         if mult == nil
 203           mult = len
 204         elsif mult > len
 205           raise "The given molecule has fewer dummy atoms (#{len}) than the bond order (#{mult})"
 206         else
 207           md.slice!(mult..-1)
 208         end
 209         sd = self.find_dummy_atoms
 210         len2 = sd.length
 211         if len2 < mult
 212           raise "Self has too few dummy atoms (#{len2}) to connect the given molecule (#{mult} bond order)"
 213         elsif len2 > mult
 214           sd.shift  #  Keep the first dummy atom untouched
 215           sd.slice!(mult..-1)
 216         end
 217         #  Find root atoms and the bond direction (which is defined by the average position of all dummy atoms)
 218         mroot = nil
 219         sroot = nil
 220         mvec = Vector3D[0, 0, 0]
 221         svec = Vector3D[0, 0, 0]
 222         md.each_index { |i|
 223           mr = mol.root_atom(md[i])
 224           if mroot == nil
 225             mroot = mr
 226           elsif mroot != mr
 227             raise "The given molecule has multiple dummy atoms bonded to different root atoms"
 228           end
 229           mvec += mol.atoms[md[i]].r
 230           sr = self.root_atom(sd[i])
 231           if sroot == nil
 232             sroot = sr
 233           elsif sroot != sr
 234             raise "self has multiple dummy atoms bonded to different root atoms"
 235           end
 236           svec += self.atoms[sd[i]].r
 237         }
 238         mvec *= 1.0 / mult
 239         svec *= 1.0 / mult
 240         mvec -= mol.atoms[mroot].r
 241         svec -= self.atoms[sroot].r
 242         #  Specify the atoms to define dihedral angles while docking
 243         #  (Preference: "untouched" dummy atoms > normal atoms > "removed" dummy atoms)
 244         mdihed = mol.atoms[mroot].connects.sort_by { |n|
 245           if md.index(n) != nil
 246             n - natoms_m
 247           elsif mol.atoms[n].element == "Du"
 248             n + natoms_m
 249           else
 250             n
 251           end
 252         } [-1]
 253         sdihed = self.atoms[sroot].connects.sort_by { |n|
 254           if sd.index(n) != nil
 255             n + natoms_s
 256           elsif self.atoms[n].element == "Du"
 257             n - natoms_s
 258           else
 259             n
 260           end
 261         } [0]
 262
 263         #  Determine the bond length
 264         me = mol.atoms[mroot].element
 265         se = self.atoms[sroot].element
 266         len = guess_bond_length(me, se, mult)
 267
 268         #  Kludge: special treatment for the amide bond
 269         if (mult == 1 && ((se == "N" && mol.is_carbonyl(mroot)) || (me == "N" && self.is_carbonyl(sroot))))
 270           self.make_sp2(sroot, sd, svec)
 271           mol.make_sp2(mroot, md, mvec)
 272           len = 1.33
 273         end
 274         if (mult == 2 && se == "C")
 275           self.make_sp2(sroot, sd, svec)
 276         end
 277         if (mult == 2 && me == "C")
 278           mol.make_sp2(mroot, md, mvec)
 279         end
 280
 281         #  Dock the two fragments
 282         self.dock(mol, sroot, mroot, svec, mvec, len, 180.0, sdihed, mdihed)
 283         #  Remove the dummy atoms
 284         md.map! { |item| item + natoms_s }
 285         self.remove(sd + md)
 286         #  Cache the list of dummy atoms
 287         @dummies = self.find_dummy_atoms
 288         self
 289   end
 290
 291   #  Internal subroutine for from_formula
 292   def Molecule.from_formula_sub(str, idx)
 293     idx0 = idx
 294     f = nil
 295         c = str[idx]
 296         if c == ?(
 297           f, idx = from_formula_sub(str, idx + 1)
 298           if str[idx] != ?)
 299             raise "Missing right parenthesis"
 300           end
 301           idx += 1
 302           if (c = str[idx]) == nil || c < ?0 || c > ?9
 303             #  No connection across the right parenthesis
 304             return f, idx
 305           end
 306         else
 307           s = str[idx..-1]
 308           key = lookup_fragment_from_string(s)
 309           if key != nil
 310             f = known_fragment(key)
 311                 idx += key.length
 312           end
 313           if f == nil
 314             ss = str[0..idx - 1] + "<?>" + str[idx..-1]
 315             raise "Cannot find atom/fragment that matches \"#{s}\": #{ss}"
 316           end
 317         end
 318         c = str[idx]
 319         valence = f.dummies.length
 320         if valence == 2 && c != nil && c >= ?0 && c <= ?9
 321           #  Concatenate the fragment n times (like (CH2)8)
 322           n = c - ?0
 323           idx += 1
 324           while (c = str[idx]) != nil && c >= ?0 && c <= ?9
 325                 n = n * 10 + c - ?0
 326                 idx += 1
 327           end
 328           if n == 0
 329                 raise "0 times repeat is not allowed"
 330           end
 331           if n > 1
 332                 f1 = f.dup
 333                 (n - 1).times { f.add_formula(f1, 1) }
 334           end
 335         end
 336         #  Add the fragment(s) from the following substring
 337         while true
 338           if (c = str[idx]) == nil || c == ?)
 339             #  End of string or parenthesized substring
 340             return f, idx
 341           end
 342           if (valence = f.dummies.length) == 0
 343                 #  This fragment is saturated
 344                 return f, idx
 345           elsif valence == 1 && idx0 > 0
 346             #  This fragment has only one dummy bond and is not the first one
 347                 return f, idx
 348           end
 349           ff, idx = from_formula_sub(str, idx)
 350           if (c = str[idx]) != nil && (c >= ?0 && c <= ?9)
 351                 #  The same fragment is added multiple times
 352             n = c - ?0
 353             idx += 1
 354             while (c = str[idx]) != nil && c >= ?0 && c <= ?9
 355                   n = n * 10 + c - ?0
 356                   idx += 1
 357             end
 358             if n == 0
 359                   raise "0 times repeat is not allowed"
 360             end
 361           else
 362             n = 1
 363           end
 364           n.times { f.add_formula(ff) }
 365         end # Repeat until the string ends or the fragment gets saturated
 366   end
 367
 368   #  Convert a string to a molecule.
 369   def Molecule.from_formula(str, resname = nil)
 370     f, idx = from_formula_sub(str, 0)
 371         f.guess_names
 372 #       f.guess_types
 373 #       if f.nresidues != 2 || resname != nil
 374           resname = "RES" if resname == nil
 375           f.assign_residue(f.all, "#{resname}.1")
 376 #       end
 377         f
 378   end
 379
 380   #  Create a fragment from str and replace the group with the generated fragment.
 381   #  This method is invoked from MainView when a detachable selection is double-clicked
 382   #  and user enters the formula in the dialog box.
 383   def dock_formula(str, group = selection)
 384     if group.length == 1
 385           #  Check if the string is an element name
 386           Parameter.atoms.each { |par|
 387                 if par.atomic_number > 0 && par.name == str
 388                   ap = atoms[group[0]]
 389                   if ap.atomic_number != par.atomic_number
 390                         if ap.name.index(ap.element) == 0
 391                           ap.name = ap.name.sub(ap.element, par.name)[0..3]
 392                         end
 393                         ap.atomic_number = par.atomic_number
 394                         ap.atom_type = ""
 395                         guess_types(group)
 396                   end
 397                   return
 398                 end
 399           }
 400         end
 401         mol = Molecule.from_formula(str)
 402         dock_fragment(mol, group)
 403   end
 404
 405   #  Replace the specified group with the given fragment.
 406   def dock_fragment(fragment, group = selection)
 407     n0 = self.natoms
 408         return if fragment.natoms == 0
 409     if group.length == 0
 410           add(fragment)
 411           sel = atom_group(n0...n0 + fragment.natoms)
 412         else
 413           frag = fragment.dup
 414       bonds_old = bonds_on_border(group)
 415           bonds_new = frag.dummies
 416           if bonds_new == nil
 417             bonds_new = frag.find_dummy_atoms()
 418           end
 419 #         puts "group = #{group.inspect}, bonds_old = #{bonds_old.inspect}, bonds_new = #{bonds_new.inspect}"
 420 #         frag.dump
 421           if bonds_old.length == 0 || bonds_new.length == 0
 422             #  Translate the added fragment to the center of the selection
 423                 v1 = self.center_of_mass(group) - frag.center_of_mass
 424                 frag.translate(v1)
 425                 remove(group)
 426                 n0 = self.natoms
 427                 add(frag)
 428                 sel = atom_group(n0...n0 + frag.natoms)
 429           else
 430             old1 = bonds_old[0][1]  #  The first "root" atom in the molecule
 431                 old2 = bonds_old[0][0]  #  The first "dummy" atom in the molecule
 432                 new1 = frag.root_atom(bonds_new[0])  #  The first "root" atom in the fragment
 433                 new2 = bonds_new[0]                  #  The first "dummy" atom in the fragment
 434                 oldv = atoms[old2].r - atoms[old1].r
 435                 newv = frag.atoms[new2].r - frag.atoms[new1].r
 436 #               puts "old1 = #{old1}, old2 = #{old2}, new1 = #{new1}, new2 = #{new2}, oldv = #{oldv.inspect}, newv = #{newv.inspect}"
 437
 438                 #  Find the most substituted atoms (which will be arranged in trans conformation)
 439                 if (atoms[old1].connects.length >= 2)
 440                   olddi1 = atoms[old1].connects.sort_by { |n| (n == old2 ? 0 : atoms[n].connects.length) } [0]
 441                 else
 442                   olddi1 = nil
 443                 end
 444                 if (frag.atoms[new1].connects.length >= 2)
 445                   newdi1 = frag.atoms[new1].connects.sort_by { |n| (n == new2 ? 0 : frag.atoms[n].connects.length) } [0]
 446                 else
 447                   newdi1 = nil
 448                 end
 449         #       dump([old1, old2, olddi1])
 450         #       frag.dump([new1, new2, newdi1])
 451                 sel = dock(frag, old1, new1, oldv, newv, 1.5, 180.0, olddi1, newdi1)  #  Add (without removing atoms)
 452                 (1..bonds_old.length - 1).each { |i|
 453                   break if i >= bonds_new.length
 454                   #  Create bonds between other "root" atoms
 455                   create_bond(bonds_old[i][1], frag.root_atom(bonds_new[i]) + n0)
 456                 }
 457                 rem1 = frag.atom_group(bonds_new)
 458                 rem = group + rem1.offset(n0)            #  The atoms to be removed
 459                 sel = atom_group(n0...n0 + frag.natoms) - rem1.offset(n0)  #  The atoms to be selected
 460                 sel = sel.deconvolute(atom_group - rem)  #  Renumber by skipping the atoms in rem
 461                 remove(rem)
 462           end
 463         end
 464         guess_types(sel)
 465         set_undoable_selection(sel)
 466   end
 467
 468
 469   #  Define molecular fragment from dump string
 470   def self.fragment_from_dump(str)
 471     f = Molecule.new.from_dump(str)
 472         f.dummies = f.find_dummy_atoms
 473         return f
 474   end
 475
 476   #  Define atom fragments
 477   elements = %w(C H Li Be B C N O F Na Mg Al Si P S Cl)
 478   valences = [4,1,1,2,3,4,3,2,1,1,2,3,4,3,2,1]
 479   degrees = {"N"=>100, "O"=>110, "P"=>95, "S"=>108}
 480   elements.each_index { |i|
 481     name = elements[i]
 482         angle = degrees[name]
 483         register_fragment(name, atom_fragment(name, valences[i], nil, angle))
 484   }
 485   #  Special fragment to represent "open valence"
 486   register_fragment("_", atom_fragment("Du", 1, nil, nil))
 487
 488   #  Define common molecular fragments
 489   formula_name = %w(CO CO2 O2C CH2 CH3 Et Pr iPr Bu tBu)
 490   formula_desc = ["CO", "C(O)O_", "OC(O)_", "CH2", "CH3", "CH2CH3", "CH2CH2CH3", "CH(CH3)2", "CH2CH2CH2CH3", "C(CH3)3"]
 491   formula_name.each_index { |i|
 492     f, idx = from_formula_sub(formula_desc[i], 0)
 493         #  Convert the "open valence" atoms to the dummy atoms
 494         f.atoms.select { |ap|
 495                 ap.element == "Du" ? ap : nil
 496         }.each_with_index { |ap, j|
 497                 ap.name = sprintf("_%d", j)
 498         }
 499         f.dummies = f.find_dummy_atoms
 500         register_fragment(formula_name[i], f)
 501   }
 502   register_fragment("Me", known_fragment("CH3"))
 503   register_fragment("Ph", fragment_from_dump(%q(
 504    0 BEN.1   C1   ca   C   -0.653   0.585  -1.068 -0.128 [1,2,10]
 505    1 BEN.1   _1   ""   Du  -1.158   1.039  -1.898  0.128 [0]
 506    2 BEN.1   C2   ca   C    0.729   0.607  -1.003 -0.128 [0,3,4]
 507    3 BEN.1   H2   ha   H    1.295   1.076  -1.783  0.128 [2]
 508    4 BEN.1   C3   ca   C    1.382   0.021   0.069 -0.128 [2,5,6]
 509    5 BEN.1   H3   ha   H    2.452   0.038   0.119  0.128 [4]
 510    6 BEN.1   C4   ca   C    0.651  -0.586   1.076 -0.128 [4,7,8]
 511    7 BEN.1   H4   ha   H    1.156  -1.039   1.906  0.128 [6]
 512    8 BEN.1   C5   ca   C   -0.732  -0.607   1.012 -0.128 [6,9,10]
 513    9 BEN.1   H5   ha   H   -1.298  -1.077   1.792  0.128 [8]
 514   10 BEN.1   C6   ca   C   -1.384  -0.021  -0.060 -0.128 [0,8,11]
 515   11 BEN.1   H6   ha   H   -2.455  -0.038  -0.110  0.128 [10])))
 516   register_fragment("C6H6", fragment_from_dump(%q(
 517    0 BEN.1   C1   ca   C   -0.653   0.585  -1.068 -0.128 [1,2,10]
 518    1 BEN.1   H1   ha   H   -1.158   1.039  -1.898  0.128 [0]
 519    2 BEN.1   C2   ca   C    0.729   0.607  -1.003 -0.128 [0,3,4]
 520    3 BEN.1   H2   ha   H    1.295   1.076  -1.783  0.128 [2]
 521    4 BEN.1   C3   ca   C    1.382   0.021   0.069 -0.128 [2,5,6]
 522    5 BEN.1   H3   ha   H    2.452   0.038   0.119  0.128 [4]
 523    6 BEN.1   C4   ca   C    0.651  -0.586   1.076 -0.128 [4,7,8]
 524    7 BEN.1   H4   ha   H    1.156  -1.039   1.906  0.128 [6]
 525    8 BEN.1   C5   ca   C   -0.732  -0.607   1.012 -0.128 [6,9,10]
 526    9 BEN.1   H5   ha   H   -1.298  -1.077   1.792  0.128 [8]
 527   10 BEN.1   C6   ca   C   -1.384  -0.021  -0.060 -0.128 [0,8,11]
 528   11 BEN.1   H6   ha   H   -2.455  -0.038  -0.110  0.128 [10])))
 529   register_fragment("C6H5", known_fragment("Ph"))
 530   register_fragment("C6H4", fragment_from_dump(%q(
 531    0 BEN.1   C1   ca   C   -0.653   0.585  -1.068 -0.128 [1,2,10]
 532    1 BEN.1   _1   ""   Du  -1.158   1.039  -1.898  0.128 [0]
 533    2 BEN.1   C2   ca   C    0.729   0.607  -1.003 -0.128 [0,3,4]
 534    3 BEN.1   H2   ha   H    1.295   1.076  -1.783  0.128 [2]
 535    4 BEN.1   C3   ca   C    1.382   0.021   0.069 -0.128 [2,5,6]
 536    5 BEN.1   H3   ha   H    2.452   0.038   0.119  0.128 [4]
 537    6 BEN.1   C4   ca   C    0.651  -0.586   1.076 -0.128 [4,7,8]
 538    7 BEN.1   _2   ""   Du   1.156  -1.039   1.906  0.128 [6]
 539    8 BEN.1   C5   ca   C   -0.732  -0.607   1.012 -0.128 [6,9,10]
 540    9 BEN.1   H5   ha   H   -1.298  -1.077   1.792  0.128 [8]
 541   10 BEN.1   C6   ca   C   -1.384  -0.021  -0.060 -0.128 [0,8,11]
 542   11 BEN.1   H6   ha   H   -2.455  -0.038  -0.110  0.128 [10])))
 543   register_fragment("C6H3", fragment_from_dump(%q(
 544    0 BEN.1   C1   ca   C   -0.653   0.585  -1.068 -0.128 [1,2,10]
 545    1 BEN.1   _1   ""   Du  -1.158   1.039  -1.898  0.128 [0]
 546    2 BEN.1   C2   ca   C    0.729   0.607  -1.003 -0.128 [0,3,4]
 547    3 BEN.1   H2   ha   H    1.295   1.076  -1.783  0.128 [2]
 548    4 BEN.1   C3   ca   C    1.382   0.021   0.069 -0.128 [2,5,6]
 549    5 BEN.1   _2   ""   Du   2.452   0.038   0.119  0.128 [4]
 550    6 BEN.1   C4   ca   C    0.651  -0.586   1.076 -0.128 [4,7,8]
 551    7 BEN.1   H4   ha   H    1.156  -1.039   1.906  0.128 [6]
 552    8 BEN.1   C5   ca   C   -0.732  -0.607   1.012 -0.128 [6,9,10]
 553    9 BEN.1   _3   ""   Du  -1.298  -1.077   1.792  0.128 [8]
 554   10 BEN.1   C6   ca   C   -1.384  -0.021  -0.060 -0.128 [0,8,11]
 555   11 BEN.1   H6   ha   H   -2.455  -0.038  -0.110  0.128 [10])))
 556
 557   #  Returns an arbitrary unit vector that is orthogonal to v
 558   def orthogonal_vector(v)
 559     vx = v.cross(Vector3D[1, 0, 0])
 560         if vx.length < 1e-3
 561           vx = v.cross(Vector3D[0, 1, 0])
 562         end
 563         vx.normalize
 564   end
 565
 566   #  Add missing hydrogen (or other atom) according to the given geometry type.
 567   #  atype can be one of the following:
 568   #  "td": tetrahedral, "tr": trigonal, "py": pyramidal (like amine nitrogen), "li": linear
 569   def add_hydrogen(idx, atype, bond = 1.07, anum = 1)
 570     nc = atoms[idx].connects.length
 571         p = atoms[idx].cr
 572         if nc == 0
 573           vx = Vector3D[1, 0, 0]
 574         else
 575           v = atoms[idx].connects.map { |i| (atoms[i].cr - p).normalize }
 576           vx = Vector3D[0, 0, 0]
 577           v.each { |v0| vx += v0 }
 578           if nc == 3
 579             if vx.length > 1e-3
 580                   vx = vx.normalize
 581                 else
 582                   #  The three atoms are in trigonal positions
 583                   vx = v[0].cross(v[1]).normalize
 584                   if vx.length < 1e-3
 585                     vx = orthogonal_vector(v[0])
 586                   end
 587                 end
 588           elsif nc == 2
 589             if vx.length < 1e-3
 590                   vx = orthognal_vector(v[0])
 591                 else
 592                   vx = vx.normalize
 593                 end
 594             vy = v[0] - v[1]
 595                 if vy.length < 1e-3
 596                   vy = orthogonal_vector(vx)
 597                 else
 598                   vy = vy.normalize
 599                 end
 600                 vz = vx.cross(vy).normalize
 601           elsif nc == 1
 602             i = atoms[idx].connects[0]
 603             #  Find most substituted atom connected to atoms[i]
 604                 vy = nil
 605                 if (atoms[i].connects.length >= 2)
 606                   j = atoms[i].connects.sort_by { |n| (n == i ? 0 : atoms[n].connects.length) }
 607                   vy = nil
 608                   j.each { |j0|
 609                     vy = vx.cross(atoms[j0].cr - atoms[i].cr)
 610                     if vy.length > 1e-3
 611                           #  The atoms j[0]-i-idx are not linear
 612                       vz = vx.cross(vy).normalize
 613                           vy = vy.normalize
 614                           break
 615                     end
 616                   }
 617                 end
 618             if !vy
 619                   vy = orthogonal_vector(vx)
 620                 else
 621                   vy = vy.normalize
 622                 end
 623                 vz = vx.cross(vy).normalize
 624           else
 625             vx = Vector3D[1, 0, 0]
 626                 vy = Vector3D[0, 1, 0]
 627                 vz = Vector3D[0, 0, 1]
 628           end
 629         end
 630         vn = []
 631         type = nil
 632     if atype == "td"
 633           raise "The atom #{idx} already has #{nc} bonds" if nc >= 4
 634           cs = -0.333333  #  cos(109.47)
 635           sn =  0.942809  #  sin(109.47) = sqrt(2)*2/3
 636           cs2 = -0.577359 #  cos(125.27)
 637           sn2 =  0.816490 #  sin(125.27)
 638           if nc == 3
 639                 vn << -vx
 640           elsif nc == 2
 641             vn << vx * cs2 + vz * sn2
 642                 vn << vx * cs2 - vz * sn2
 643           else
 644             vn << vx if nc == 0
 645             vn << vx * cs + vy * sn
 646                 vn << vx * cs - vy * sn * 0.5 + vz * sn * 0.86603
 647                 vn << vx * cs - vy * sn * 0.5 - vz * sn * 0.86603
 648           end
 649           type = "hc" if anum == 1
 650         elsif atype == "tr"
 651           raise "The atom #{idx} already has #{nc} bonds" if nc >= 3
 652           if nc == 2
 653             vn << -vx
 654           else
 655             vn << vx if nc == 0
 656                 vn << vx * 0.5 + vy * 0.86603
 657                 vn << vx * 0.5 - vy * 0.86603
 658           end
 659           type = "ha" if anum == 1
 660         elsif atype == "py"
 661           raise "The atom #{idx} already has #{nc} bonds" if nc >= 3
 662           cs = -0.292376  #  cos(107)
 663           sn =  0.956303  #  sin(107)
 664           cs2 = -0.491527 #  cos(119.4)
 665           sn2 =  0.870862 #  sin(119.4)
 666           if nc == 2
 667             vn << vx * cs2 + vz * sn2
 668           else
 669             vn << vx if nc == 0
 670                 vn << vx * cs + vy * sn
 671                 vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt(1 + 2 * cs))
 672           end
 673           type = "h2" if anum == 1
 674         elsif atype == "be"
 675           raise "The atom #{idx} already has #{nc} bonds" if nc >= 2
 676           cs = -0.258819  #  cos(105)
 677           sn =  0.965926  #  sin(105)
 678           vn << vx if nc == 0
 679           vn << vx * cs + vy * sn
 680           type = "ho" if anum == 1
 681         elsif atype == "li"
 682           raise "The atom #{idx} already has #{nc} bonds" if nc >= 2
 683           vn << vx if nc == 0
 684           vn << -vx
 685           type = "hz" if anum == 1
 686         else
 687           raise "Unknown atom type #{atype}"
 688         end
 689         aname = Parameter.atom(anum).name
 690         name = atoms[idx].name
 691         name = name.gsub(/\A[A-Za-z]*/, aname)[0..2]
 692         vn.each_with_index { |v, i|
 693           ap = create_atom(sprintf("%s%d", name, i + 1), idx + i + 1)
 694           ap.atomic_number = anum
 695           ap.atom_type = (type || ap.element)
 696           assign_residue(atom_group(ap.index), atoms[idx].res_seq)
 697           ap.r = atoms[idx].r + v * bond
 698           create_bond(idx, ap.index)
 699         }
 700         #  Update selection
 701 #       if selection.count > 0
 702 #         sel = selection.convolute(IntGroup[0..idx, (idx + vn.count + 1)..(natoms - 1)])
 703 #         self.selection = sel
 704 #       end
 705   end
 706
 707   def add_hydrogen_on_group(group, atype, bond = 1.07, anum = 1)
 708     group.reverse_each { |i|
 709           add_hydrogen(i, atype, bond, anum)
 710         }
 711     self
 712   end
 713
 714 end