4 # Created by Toshi Nagata.
5 # Copyright 2008 Toshi Nagata. All rights reserved.
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.
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.
18 @@known_fragments = Hash.new
19 @@known_fragment_names = []
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
28 attr_accessor :dummies
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)
35 if ap.element == "Du" && ap.name[0] == ?_
37 if how_many != nil && d.length >= how_many
45 # Returns the first atom connected to the specified atom.
46 def root_atom(dummy_atom)
47 return atoms[dummy_atom].connects[0]
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
60 @@known_fragment_names[idx, 0] = name
63 def Molecule.known_fragment(name)
64 @@known_fragments[name].dup
67 def Molecule.lookup_fragment_from_string(str)
68 @@known_fragment_names.each { |name|
69 if str.rindex(name, 0) == 0
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
84 fragment.add_atom("_1", "", "Du", length, 0)
87 degree = 180.0 if valence == 2
88 degree = 120.0 if valence == 3
89 degree = 109.2 if valence == 4
91 fragment.add_atom("_2", "", "Du", length, 0, degree, 1)
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)
100 dihed = Math.acos(acs) * Rad2Deg
105 fragment.add_atom("_3", "", "Du", length, 0, degree, 1, dihed, 2)
107 fragment.add_atom("_4", "", "Du", length, 0, degree, 1, -dihed, 2)
112 fragment.dummies = fragment.find_dummy_atoms
116 # Examine whether the given atom is a carbonyl carbon
118 if atoms[idx].element == "C"
119 atoms[idx].connects.each { |idx2|
120 if atoms[idx2].element == "O" && atoms[idx2].connects.length == 1
128 def guess_bond_length(e1, e2, mult)
134 if (e1 == "C" && e2 == "C")
135 len = (mult == 3 ? 1.18 : (mult == 2 ? 1.39 : 1.54))
136 elsif (e1 == "C" && e2 == "H")
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")
144 elsif (e1 == "H" && e2 == "N")
146 elsif (e1 == "C" && e2 == "Cl")
148 elsif (e1 == "C" && e2 == "S")
150 elsif (e1 == "O" && e2 == "S")
151 len = (mult == 2 ? 1.60 : 1.80)
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
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")
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)
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
200 raise "The given molecule has no dummy atoms"
205 raise "The given molecule has fewer dummy atoms (#{len}) than the bond order (#{mult})"
209 sd = self.find_dummy_atoms
212 raise "Self has too few dummy atoms (#{len2}) to connect the given molecule (#{mult} bond order)"
214 sd.shift # Keep the first dummy atom untouched
217 # Find root atoms and the bond direction (which is defined by the average position of all dummy atoms)
220 mvec = Vector3D[0, 0, 0]
221 svec = Vector3D[0, 0, 0]
223 mr = mol.root_atom(md[i])
227 raise "The given molecule has multiple dummy atoms bonded to different root atoms"
229 mvec += mol.atoms[md[i]].r
230 sr = self.root_atom(sd[i])
234 raise "self has multiple dummy atoms bonded to different root atoms"
236 svec += self.atoms[sd[i]].r
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
247 elsif mol.atoms[n].element == "Du"
253 sdihed = self.atoms[sroot].connects.sort_by { |n|
254 if sd.index(n) != nil
256 elsif self.atoms[n].element == "Du"
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)
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)
274 if (mult == 2 && se == "C")
275 self.make_sp2(sroot, sd, svec)
277 if (mult == 2 && me == "C")
278 mol.make_sp2(mroot, md, mvec)
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 }
286 # Cache the list of dummy atoms
287 @dummies = self.find_dummy_atoms
291 # Internal subroutine for from_formula
292 def Molecule.from_formula_sub(str, idx)
297 f, idx = from_formula_sub(str, idx + 1)
299 raise "Missing right parenthesis"
302 if (c = str[idx]) == nil || c < ?0 || c > ?9
303 # No connection across the right parenthesis
308 key = lookup_fragment_from_string(s)
310 f = known_fragment(key)
314 ss = str[0..idx - 1] + "<?>" + str[idx..-1]
315 raise "Cannot find atom/fragment that matches \"#{s}\": #{ss}"
319 valence = f.dummies.length
320 if valence == 2 && c != nil && c >= ?0 && c <= ?9
321 # Concatenate the fragment n times (like (CH2)8)
324 while (c = str[idx]) != nil && c >= ?0 && c <= ?9
329 raise "0 times repeat is not allowed"
333 (n - 1).times { f.add_formula(f1, 1) }
336 # Add the fragment(s) from the following substring
338 if (c = str[idx]) == nil || c == ?)
339 # End of string or parenthesized substring
342 if (valence = f.dummies.length) == 0
343 # This fragment is saturated
345 elsif valence == 1 && idx0 > 0
346 # This fragment has only one dummy bond and is not the first one
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
354 while (c = str[idx]) != nil && c >= ?0 && c <= ?9
359 raise "0 times repeat is not allowed"
364 n.times { f.add_formula(ff) }
365 end # Repeat until the string ends or the fragment gets saturated
368 # Convert a string to a molecule.
369 def Molecule.from_formula(str, resname = nil)
370 f, idx = from_formula_sub(str, 0)
373 # if f.nresidues != 2 || resname != nil
374 resname = "RES" if resname == nil
375 f.assign_residue(f.all, "#{resname}.1")
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)
385 # Check if the string is an element name
386 Parameter.builtin.atoms.each { |par|
387 if par.atomic_number > 0 && par.name == str
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]
393 ap.atomic_number = par.atomic_number
401 mol = Molecule.from_formula(str)
402 dock_fragment(mol, group)
405 # Replace the specified group with the given fragment.
406 def dock_fragment(fragment, group = selection)
408 return if fragment.natoms == 0
411 sel = atom_group(n0...n0 + fragment.natoms)
414 bonds_old = bonds_on_border(group)
415 bonds_new = frag.dummies
417 bonds_new = frag.find_dummy_atoms()
419 # puts "group = #{group.inspect}, bonds_old = #{bonds_old.inspect}, bonds_new = #{bonds_new.inspect}"
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
428 sel = atom_group(n0...n0 + frag.natoms)
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}"
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]
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]
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)
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
465 set_undoable_selection(sel)
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
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|
482 angle = degrees[name]
483 register_fragment(name, atom_fragment(name, valences[i], nil, angle))
485 # Special fragment to represent "open valence"
486 register_fragment("_", atom_fragment("Du", 1, nil, nil))
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)
499 f.dummies = f.find_dummy_atoms
500 register_fragment(formula_name[i], f)
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])))
557 # Returns an arbitrary unit vector that is orthogonal to v
558 def orthogonal_vector(v)
559 vx = v.cross(Vector3D[1, 0, 0])
561 vx = v.cross(Vector3D[0, 1, 0])
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
573 vx = Vector3D[1, 0, 0]
575 v = atoms[idx].connects.map { |i| (atoms[i].cr - p).normalize }
576 vx = Vector3D[0, 0, 0]
577 v.each { |v0| vx += v0 }
582 # The three atoms are in trigonal positions
583 vx = v[0].cross(v[1]).normalize
585 vx = orthogonal_vector(v[0])
590 vx = orthognal_vector(v[0])
596 vy = orthogonal_vector(vx)
600 vz = vx.cross(vy).normalize
602 i = atoms[idx].connects[0]
603 # Find most substituted atom connected to atoms[i]
605 if (atoms[i].connects.length >= 2)
606 j = atoms[i].connects.sort_by { |n| (n == i ? 0 : atoms[n].connects.length) }
609 vy = vx.cross(atoms[j0].cr - atoms[i].cr)
611 # The atoms j[0]-i-idx are not linear
612 vz = vx.cross(vy).normalize
619 vy = orthogonal_vector(vx)
623 vz = vx.cross(vy).normalize
625 vx = Vector3D[1, 0, 0]
626 vy = Vector3D[0, 1, 0]
627 vz = Vector3D[0, 0, 1]
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)
641 vn << vx * cs2 + vz * sn2
642 vn << vx * cs2 - vz * sn2
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
649 type = "hc" if anum == 1
651 raise "The atom #{idx} already has #{nc} bonds" if nc >= 3
656 vn << vx * 0.5 + vy * 0.86603
657 vn << vx * 0.5 - vy * 0.86603
659 type = "ha" if anum == 1
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)
667 vn << vx * cs2 + vz * sn2
670 vn << vx * cs + vy * sn
671 vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt(1 + 2 * cs))
673 type = "h2" if anum == 1
675 raise "The atom #{idx} already has #{nc} bonds" if nc >= 2
676 cs = -0.258819 # cos(105)
677 sn = 0.965926 # sin(105)
679 vn << vx * cs + vy * sn
680 type = "ho" if anum == 1
682 raise "The atom #{idx} already has #{nc} bonds" if nc >= 2
685 type = "hz" if anum == 1
687 raise "Unknown atom type #{atype}"
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)
701 # if selection.count > 0
702 # sel = selection.convolute(IntGroup[0..idx, (idx + vn.count + 1)..(natoms - 1)])
703 # self.selection = sel
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)