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

diff --git a/Scripts/crystal.rb b/Scripts/crystal.rb
index 2a94ef8..4d0d748 100755 (executable)
--- a/Scripts/crystal.rb
+++ b/Scripts/crystal.rb
@@ -29,19 +29,6 @@ class AtomRef
   end
 end
 
-#  Used in bond_angle_with_sigma
-module Math
-  def acos_safe(arg)
-    if arg <= -1.0
-      return PI
-    elsif arg >= 1.0
-      return 0.0
-    else
-      return acos(arg)
-    end
-  end
-end
-
 module Kernel
 
 def symmetry_to_string(sym)
@@ -436,12 +423,6 @@ end
 
 end
 
-module Math
-  def sqrt_safe(arg)
-    arg <= 0.0 ? 0.0 : sqrt(arg)
-  end
-end
-
 #  Best-fit planes
 #  Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
 #       T. Ito, Acta Cryst 1981, A37, 621-624
@@ -1826,7 +1807,7 @@ end
                bv = Vector3D[cos(gamma), sin(gamma), 0]
                cx = cos(beta)
                cy = (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma)
-               cz = sqrt(1.0 - cx * cx - cy * cy)
+               cz = sqrt_safe(1.0 - cx * cx - cy * cy)
                cv = Vector3D[cx, cy, cz]
                x0 = @box[0].normalize
                z0 = (@box[0].cross(@box[1])).normalize

diff --git a/Scripts/formula.rb b/Scripts/formula.rb
index c44d650..6e2f496 100755 (executable)
--- a/Scripts/formula.rb
+++ b/Scripts/formula.rb
@@ -641,7 +641,7 @@ class Molecule
          else
            vn << vx if nc == 0
                vn << vx * cs + vy * sn
-               vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt(1 + 2 * cs))
+               vn << vx * cs + vy * (cs * (1 - cs) / sn) + vz * ((1 - cs) / sn * Math.sqrt_safe(1 + 2 * cs))
          end
          type = "h2" if anum == 1
        elsif atype == "be"

diff --git a/Scripts/md.rb b/Scripts/md.rb
index b68c111..edd0766 100755 (executable)
--- a/Scripts/md.rb
+++ b/Scripts/md.rb
@@ -1038,7 +1038,7 @@ class Molecule
           eps = pp.eps
           d = pp.r_eq * 2
         else
-          eps = Math.sqrt(p1.eps * p2.eps)
+          eps = Math.sqrt_safe(p1.eps * p2.eps)
           d = p1.r_eq + p2.r_eq
         end
         vdw_access_table[i * nvdw_pars + j] =  vdw_access_table[j * nvdw_pars + i] = vdw_pair_table.length

diff --git a/Scripts/molecule.rb b/Scripts/molecule.rb
index 0d1817e..b8780e7 100755 (executable)
--- a/Scripts/molecule.rb
+++ b/Scripts/molecule.rb
@@ -129,7 +129,7 @@ class Molecule
       return 0.0
     end
     cs = r21.dot(r23) / (w1 * w2)
-    Math.atan2(Math.sqrt(1 - cs*cs), cs) * Rad2Deg
+    Math.atan2(Math.sqrt_safe(1 - cs*cs), cs) * Rad2Deg
   end
   
   #  Calculate the dihedral angle defined by four vectors.
@@ -280,7 +280,10 @@ class Molecule
       else
         v3 = v2.cross(v1).normalize
       end
-      angle = Math.atan2(Math.sqrt(1.0 - cs*cs), cs) * Rad2Deg
+         #if cs > 1.0
+           puts "cs = #{cs}"
+         #end
+      angle = Math.atan2(Math.sqrt_safe(1.0 - cs*cs), cs) * Rad2Deg
       mol.rotate(v3, angle, mol.atoms[base2].r)
     end
     #  Move the second molecule so that the atom 'base2' is located at atoms[base1].r+v1*len

diff --git a/Scripts/startup.rb b/Scripts/startup.rb
index 0a55a35..06f7fb7 100755 (executable)
--- a/Scripts/startup.rb
+++ b/Scripts/startup.rb
@@ -57,6 +57,23 @@ module Enumerable
   end
 end
 
+module Math
+  def acos_safe(arg)
+    if arg <= -1.0
+      return PI
+    elsif arg >= 1.0
+      return 0.0
+    else
+      return acos(arg)
+    end
+  end
+
+  def sqrt_safe(arg)
+    arg <= 0.0 ? 0.0 : sqrt(arg)
+  end
+end
+
+
 module Kernel
   def filecopy(src, dst)
     fpin = File.open(src, "rb")

diff --git a/Scripts/transform.rb b/Scripts/transform.rb
index d1adceb..9ab45aa 100755 (executable)
--- a/Scripts/transform.rb
+++ b/Scripts/transform.rb
@@ -34,25 +34,25 @@ class Transform
     zz = self[2,2]
     ww = xx + yy + zz + 1.0
     if (ww >= 1.0)
-      w = Math.sqrt(0.25 * ww)
+      w = Math.sqrt_safe(0.25 * ww)
       ww = 0.25 / w;
       x = (self[1,2] - self[2,1]) * ww
       y = (self[2,0] - self[0,2]) * ww
       z = (self[0,1] - self[1,0]) * ww
     elsif (xx >= yy && xx >= zz)
-      x = Math.sqrt(0.25 * (1.0 + xx - yy - zz))
+      x = Math.sqrt_safe(0.25 * (1.0 + xx - yy - zz))
       ww = 0.25 / x
       y = (self[1,0] + self[0,1]) * ww
       z = (self[2,0] + self[0,2]) * ww
       w = (self[1,2] - self[2,1]) * ww
     elsif (yy >= zz)
-      y = Math.sqrt(0.25 * (1.0 + yy - xx - zz))
+      y = Math.sqrt_safe(0.25 * (1.0 + yy - xx - zz))
       ww = 0.25 / y
       x = (self[1,0] + self[0,1]) * ww
       z = (self[1,2] + self[2,1]) * ww
       w = (self[2,0] - self[0,2]) * ww
     else
-      z = Math.sqrt(0.25 * (1.0 + zz - xx - yy))
+      z = Math.sqrt_safe(0.25 * (1.0 + zz - xx - yy))
       ww = 0.25 / z
       x = (self[2,0] + self[0,2]) * ww
       y = (self[2,1] + self[1,2]) * ww

diff --git a/Scripts/uff.rb b/Scripts/uff.rb
index 0528585..e0477fa 100644 (file)
--- a/Scripts/uff.rb
+++ b/Scripts/uff.rb
@@ -140,8 +140,8 @@ UFFParams = [
 #  Calculate UFF bond length
 def uff_bond_length(r1, r2, x1, x2, bond_order)
   bond_order_correction = -0.1332 * (r1 + r2) * log(bond_order)
-  sq = sqrt(x1) - sqrt(x2)
-  electronegativity_correction = r1 * r2 * (sqrt(x1) - sqrt(x2)) ** 2 / (x1 * r1 + x2 * r2)
+  sq = sqrt_safe(x1) - sqrt_safe(x2)
+  electronegativity_correction = r1 * r2 * (sqrt_safe(x1) - sqrt_safe(x2)) ** 2 / (x1 * r1 + x2 * r2)
   return r1 + r2 + bond_order_correction + electronegativity_correction
 end
 
@@ -161,7 +161,7 @@ def uff_angle_force(idx1, idx2, idx3, bond_order_12, bond_order_23, angle)
   cost = cos(angle * 3.1415927 / 180.0)
   r12 = uff_bond_length(r1, r2, UFFParams[idx1][10], UFFParams[idx2][10], bond_order_12)
   r23 = uff_bond_length(r2, r3, UFFParams[idx2][10], UFFParams[idx3][10], bond_order_23)
-  r13 = sqrt(r12 * r12 + r23 * r23 - 2 * r12 * r23 * cost)
+  r13 = sqrt_safe(r12 * r12 + r23 * r23 - 2 * r12 * r23 * cost)
   return 332.06 * UFFParams[idx1][7] * UFFParams[idx3][7] / (r13 ** 5) * (r12 * r23 * (1.0 - cost * cost) - r13 * r13 * cost)
 end