OSDN -- オープンソース・ソフトウェアの開発とダウンロード
日本語の翻訳状況 translation status for ja
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Molby: コミット

Molecular Modeling Software


コミットメタ情報

リビジョン8ec1148d50b4484e30ae44b4ba70813fa60522dd (tree)
日時2014-04-03 23:41:54
作者toshinagata1964 <toshinagata1964@a2be...>
コミッターtoshinagata1964

ログメッセージ

sqrt out-of-range error is now avoided in docking a fragment

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@524 a2be9bc6-48de-4e38-9406-05402d4bc13c

変更サマリ

差分

--- a/Scripts/crystal.rb
+++ b/Scripts/crystal.rb
@@ -29,19 +29,6 @@ class AtomRef
2929 end
3030 end
3131
32-# Used in bond_angle_with_sigma
33-module Math
34- def acos_safe(arg)
35- if arg <= -1.0
36- return PI
37- elsif arg >= 1.0
38- return 0.0
39- else
40- return acos(arg)
41- end
42- end
43-end
44-
4532 module Kernel
4633
4734 def symmetry_to_string(sym)
@@ -436,12 +423,6 @@ end
436423
437424 end
438425
439-module Math
440- def sqrt_safe(arg)
441- arg <= 0.0 ? 0.0 : sqrt(arg)
442- end
443-end
444-
445426 # Best-fit planes
446427 # Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
447428 # T. Ito, Acta Cryst 1981, A37, 621-624
@@ -1826,7 +1807,7 @@ end
18261807 bv = Vector3D[cos(gamma), sin(gamma), 0]
18271808 cx = cos(beta)
18281809 cy = (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma)
1829- cz = sqrt(1.0 - cx * cx - cy * cy)
1810+ cz = sqrt_safe(1.0 - cx * cx - cy * cy)
18301811 cv = Vector3D[cx, cy, cz]
18311812 x0 = @box[0].normalize
18321813 z0 = (@box[0].cross(@box[1])).normalize
--- a/Scripts/formula.rb
+++ b/Scripts/formula.rb
@@ -641,7 +641,7 @@ class Molecule
641641 else
642642 vn << vx if nc == 0
643643 vn << vx * cs + vy * sn
644- vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt(1 + 2 * cs))
644+ vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt_safe(1 + 2 * cs))
645645 end
646646 type = "h2" if anum == 1
647647 elsif atype == "be"
--- a/Scripts/md.rb
+++ b/Scripts/md.rb
@@ -1038,7 +1038,7 @@ class Molecule
10381038 eps = pp.eps
10391039 d = pp.r_eq * 2
10401040 else
1041- eps = Math.sqrt(p1.eps * p2.eps)
1041+ eps = Math.sqrt_safe(p1.eps * p2.eps)
10421042 d = p1.r_eq + p2.r_eq
10431043 end
10441044 vdw_access_table[i * nvdw_pars + j] = vdw_access_table[j * nvdw_pars + i] = vdw_pair_table.length
--- a/Scripts/molecule.rb
+++ b/Scripts/molecule.rb
@@ -129,7 +129,7 @@ class Molecule
129129 return 0.0
130130 end
131131 cs = r21.dot(r23) / (w1 * w2)
132- Math.atan2(Math.sqrt(1 - cs*cs), cs) * Rad2Deg
132+ Math.atan2(Math.sqrt_safe(1 - cs*cs), cs) * Rad2Deg
133133 end
134134
135135 # Calculate the dihedral angle defined by four vectors.
@@ -280,7 +280,10 @@ class Molecule
280280 else
281281 v3 = v2.cross(v1).normalize
282282 end
283- angle = Math.atan2(Math.sqrt(1.0 - cs*cs), cs) * Rad2Deg
283+ #if cs > 1.0
284+ puts "cs = #{cs}"
285+ #end
286+ angle = Math.atan2(Math.sqrt_safe(1.0 - cs*cs), cs) * Rad2Deg
284287 mol.rotate(v3, angle, mol.atoms[base2].r)
285288 end
286289 # Move the second molecule so that the atom 'base2' is located at atoms[base1].r+v1*len
--- a/Scripts/startup.rb
+++ b/Scripts/startup.rb
@@ -57,6 +57,23 @@ module Enumerable
5757 end
5858 end
5959
60+module Math
61+ def acos_safe(arg)
62+ if arg <= -1.0
63+ return PI
64+ elsif arg >= 1.0
65+ return 0.0
66+ else
67+ return acos(arg)
68+ end
69+ end
70+
71+ def sqrt_safe(arg)
72+ arg <= 0.0 ? 0.0 : sqrt(arg)
73+ end
74+end
75+
76+
6077 module Kernel
6178 def filecopy(src, dst)
6279 fpin = File.open(src, "rb")
--- a/Scripts/transform.rb
+++ b/Scripts/transform.rb
@@ -34,25 +34,25 @@ class Transform
3434 zz = self[2,2]
3535 ww = xx + yy + zz + 1.0
3636 if (ww >= 1.0)
37- w = Math.sqrt(0.25 * ww)
37+ w = Math.sqrt_safe(0.25 * ww)
3838 ww = 0.25 / w;
3939 x = (self[1,2] - self[2,1]) * ww
4040 y = (self[2,0] - self[0,2]) * ww
4141 z = (self[0,1] - self[1,0]) * ww
4242 elsif (xx >= yy && xx >= zz)
43- x = Math.sqrt(0.25 * (1.0 + xx - yy - zz))
43+ x = Math.sqrt_safe(0.25 * (1.0 + xx - yy - zz))
4444 ww = 0.25 / x
4545 y = (self[1,0] + self[0,1]) * ww
4646 z = (self[2,0] + self[0,2]) * ww
4747 w = (self[1,2] - self[2,1]) * ww
4848 elsif (yy >= zz)
49- y = Math.sqrt(0.25 * (1.0 + yy - xx - zz))
49+ y = Math.sqrt_safe(0.25 * (1.0 + yy - xx - zz))
5050 ww = 0.25 / y
5151 x = (self[1,0] + self[0,1]) * ww
5252 z = (self[1,2] + self[2,1]) * ww
5353 w = (self[2,0] - self[0,2]) * ww
5454 else
55- z = Math.sqrt(0.25 * (1.0 + zz - xx - yy))
55+ z = Math.sqrt_safe(0.25 * (1.0 + zz - xx - yy))
5656 ww = 0.25 / z
5757 x = (self[2,0] + self[0,2]) * ww
5858 y = (self[2,1] + self[1,2]) * ww
--- a/Scripts/uff.rb
+++ b/Scripts/uff.rb
@@ -140,8 +140,8 @@ UFFParams = [
140140 # Calculate UFF bond length
141141 def uff_bond_length(r1, r2, x1, x2, bond_order)
142142 bond_order_correction = -0.1332 * (r1 + r2) * log(bond_order)
143- sq = sqrt(x1) - sqrt(x2)
144- electronegativity_correction = r1 * r2 * (sqrt(x1) - sqrt(x2)) ** 2 / (x1 * r1 + x2 * r2)
143+ sq = sqrt_safe(x1) - sqrt_safe(x2)
144+ electronegativity_correction = r1 * r2 * (sqrt_safe(x1) - sqrt_safe(x2)) ** 2 / (x1 * r1 + x2 * r2)
145145 return r1 + r2 + bond_order_correction + electronegativity_correction
146146 end
147147
@@ -161,7 +161,7 @@ def uff_angle_force(idx1, idx2, idx3, bond_order_12, bond_order_23, angle)
161161 cost = cos(angle * 3.1415927 / 180.0)
162162 r12 = uff_bond_length(r1, r2, UFFParams[idx1][10], UFFParams[idx2][10], bond_order_12)
163163 r23 = uff_bond_length(r2, r3, UFFParams[idx2][10], UFFParams[idx3][10], bond_order_23)
164- r13 = sqrt(r12 * r12 + r23 * r23 - 2 * r12 * r23 * cost)
164+ r13 = sqrt_safe(r12 * r12 + r23 * r23 - 2 * r12 * r23 * cost)
165165 return 332.06 * UFFParams[idx1][7] * UFFParams[idx3][7] / (r13 ** 5) * (r12 * r23 * (1.0 - cost * cost) - r13 * r13 * cost)
166166 end
167167
旧リポジトリブラウザで表示