Scripts/crystal.rb

   1 #
   2 #  crystal.rb
   3 #
   4 #  Created by Toshi Nagata.
   5 #  Copyright 2012 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 #  Definition for use in ORTEP
  17 class AtomRef
  18   def to_adc
  19     sym = self.symop
  20     if sym == nil
  21       idx = self.index + 1
  22       symcode = 55501
  23     else
  24       idx = sym[4] + 1
  25       symcode = (sym[1] + 5) * 10000 + (sym[2] + 5) * 1000 + (sym[3] + 5) * 100 + sym[0] + 1
  26     end
  27     return idx, symcode
  28   end
  29 end
  30
  31 #  Used in bond_angle_with_sigma
  32 module Math
  33   def acos_safe(arg)
  34     if arg <= -1.0
  35       return PI
  36     elsif arg >= 1.0
  37       return 0.0
  38     else
  39       return acos(arg)
  40     end
  41   end
  42 end
  43
  44 class Molecule
  45
  46 #  Export ortep to File fp
  47 #  If attr is a hash, it represents options for drawing.
  48 #  "atoms"=>[[group, type, color, rad], [group, type, color, rad], ...]
  49 #    group is an IntGroup, representing a group of atoms.
  50 #    type is an atom type: 0 = boundary, 1 = boundary + principal, 2 = 1 + axes, 3 = 2 + shades, 4 = fixed size
  51 #    color is 0..7 (0,1=black, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 7=yellow)
  52 #    rad is the radius of the atom (in Angstrom) whose type is fixed size
  53 #    If an atom appears in multiple entries, the later entry is valid.
  54 #  "bonds"=>[[group1, group2, type, color], [group1, group2, type, color], ...]
  55 #    group1 and group2 are IntGroups, reprensenting the atoms constituting the bonds.
  56 #    type is a bond type: 0 to 4 for bonds with no shades to 4 shades.
  57 #    color is 0..7 as in atoms.
  58 #    If a bond appears in multiple entries, the later entry is valid.
  59 def export_ortep(fp, attr = nil)
  60
  61   #  Create atom list
  62   hidden = atom_group { |ap| !is_atom_visible(ap.index) }
  63   hydrogen = self.show_hydrogens
  64   expanded = self.show_expanded
  65   atomlist = atom_group { |ap|
  66     (ap.element != "H" || hydrogen) &&
  67     (ap.symop == nil || expanded) &&
  68     (!hidden.include?(ap.index))
  69   }
  70
  71   #  Title
  72   fp.printf "%-78.78s\n", self.name + ": generated by Molby at " + Time.now.to_s
  73
  74   #  Cell parameters
  75   cp = self.cell
  76   if cp == nil
  77     cp = [1, 1, 1, 90, 90, 90]
  78   end
  79   fp.printf "%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f\n", cp[0], cp[1], cp[2], cp[3], cp[4], cp[5]
  80
  81   #  Symmetry operations
  82   syms = self.symmetries
  83   if syms == nil || syms.length == 0
  84    fp.print "1             0  1  0  0              0  0  1  0              0  0  0  1\n"
  85   else
  86     syms.each_with_index { |s, i|
  87       a = s.to_a
  88       fp.printf "%s%14g%3g%3g%3g%15g%3g%3g%3g%15g%3g%3g%3g\n", (i == syms.length - 1 ? "1" : " "), a[9], a[0], a[1], a[2], a[10], a[3], a[4], a[5], a[11], a[6], a[7], a[8]
  89     }
  90   end
  91
  92   #  Atoms (all symmetry unique atoms regardless they are visible or not)
  93   n = 0
  94   aattr = (attr ? attr["atoms"] : nil)
  95   each_atom { |ap|
  96     break if ap.symop != nil
  97     fp.printf " %4.4s%22s%9.4f%9.4f%9.4f%9d\n", ap.name, "", ap.fract_x, ap.fract_y, ap.fract_z, 0
  98         rad = -1.0
  99         if aattr
 100           aattr.reverse_each { |at|
 101             if at[0].include?(ap.index)
 102                   if at[1] == 4
 103                     rad = at[3].to_f
 104                     if rad == 0.0
 105                       rad = 0.1
 106                     end
 107                   end
 108                   break
 109                 end
 110           }
 111         end
 112     an = ap.aniso
 113     if an != nil && rad < 0.0
 114       fp.printf " %8.5f%9.6f%9.6f%9.6f%9.6f%9.6f%9d\n", an[0], an[1], an[2], an[3], an[4], an[5], 0
 115     else
 116           type = 7
 117           if rad < 0.0
 118         rad = ap.temp_factor
 119         rad = 1.2 if rad <= 0
 120             type = 6
 121           end
 122       fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", rad, 0.0, 0.0, 0.0, 0.0, 0.0, type
 123     end
 124     n += 1
 125   }
 126   natoms_tep = n
 127
 128   #  Special points to specify cartesian axes
 129   axis, angle = self.get_view_rotation
 130   tr = Transform.rotation(axis, -angle)
 131   org = self.get_view_center
 132   x = org + tr.column(0)
 133   y = org + tr.column(1)
 134   tr = self.cell_transform
 135   if tr
 136     tr = tr.inverse
 137   else
 138     tr = Transform.identity
 139   end
 140   org = tr * org
 141   x = tr * x
 142   y = tr * y
 143   fp.printf " CNTR                      %9.4f%9.4f%9.4f        0\n", org.x, org.y, org.z
 144   fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
 145   fp.printf " X                         %9.4f%9.4f%9.4f        0\n", x.x, x.y, x.z
 146   fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
 147   fp.printf " Y                         %9.4f%9.4f%9.4f        0\n", y.x, y.y, y.z
 148   fp.printf "1%8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
 149
 150   #  Initialize
 151   fp.print  "      201\n"
 152
 153   #  Pen size
 154   fp.print  "      205        8\n"
 155
 156   #  Paper size, margin, viewing distance
 157   fp.print  "      301      6.6      6.6        0      0.0\n"
 158
 159   #  Sort atoms by symop and index
 160   acodes = Hash.new
 161   ig = IntGroup.new   #  Used for collecting contiguous atoms
 162   each_atom(atomlist) { |ap|
 163     idx, symcode = ap.to_adc
 164     acode = symcode * self.natoms + idx - 1
 165     acodes[acode] = ap.index
 166     ig.add(acode)
 167   }
 168   index2an = []       #  Index in Molecule to Index in ORTEP
 169   ig.each_with_index { |acode, i|
 170     index2an[acodes[acode]] = i + 1
 171   }
 172   i = 0
 173   adcs = []
 174   while (r = ig.range_at(i)) != nil
 175     s = r.first
 176     s = (s / self.natoms) + (s % self.natoms + 1) * 100000  #  Rebuild true ADC (atom index is in the upper digit)
 177     e = r.last
 178     e = (e / self.natoms) + (e % self.natoms + 1) * 100000
 179     if s < e
 180       adcs.push(s)
 181       adcs.push(-e)
 182     else
 183       adcs.push(s)
 184     end
 185     i += 1
 186   end
 187   k = 0
 188
 189   #  Atom list
 190   adcs.each_with_index { |a, i|
 191     if k == 0
 192       fp.print "      401"
 193     end
 194     fp.printf "%9d", a
 195     k += 1
 196     if i == adcs.length - 1 || k == 6 || (k == 5 && i < adcs.length - 2 && adcs[i + 2] < 0)
 197       fp.print "\n"
 198       k = 0
 199     end
 200   }
 201
 202   #  Axes
 203   fp.printf "      501%4d55501%4d55501%4d55501%4d55501%4d55501                 1\n", natoms_tep + 1, natoms_tep + 1, natoms_tep + 2, natoms_tep + 1, natoms_tep + 3
 204 #  fp.print  "      502        1      0.0        2      0.0        3      0.0\n"
 205
 206   #  Autoscale
 207   fp.print  "      604                               1.538\n"
 208
 209   #  Explicit bonds
 210   bond_inst = Array.new(35) { [] }   #  Bonds for types 0 to 4 and colors 1 to 7
 211   battr = (attr ? attr["bonds"] : nil)
 212   bonds.each { |b|
 213     next if !atomlist.include?(b[0]) || !atomlist.include?(b[1])
 214         btype = bcol = nil
 215         if battr
 216           battr.reverse_each { |at|
 217             if (at[0].include?(b[0]) && at[1].include?(b[1])) || (at[0].include?(b[1]) && at[1].include?(b[0]))
 218                   btype = at[2] % 5
 219                   bcol = (at[3] != 0 ? at[3] - 1 : at[3]) % 7
 220                   break
 221                 end
 222           }
 223         end
 224         if btype == nil
 225           bcol = 0
 226       an1 = atoms[b[0]].atomic_number
 227       an2 = atoms[b[1]].atomic_number
 228       if an1 == 1 || an2 == 1
 229         btype = 0
 230       elsif an1 <= 8 && an2 <= 8
 231         btype = 2
 232       else
 233         btype = 4
 234       end
 235         end
 236     bond_inst[btype * 7 + bcol].push(b[0], b[1])
 237   }
 238
 239   #  Output bond specifications
 240   #  Avoid including too many ADCs in a single 811/821 instruction
 241   #  (Upper limit is 140. Here we divide at every 36 ADCs)
 242   output_bonds = lambda { |icode|
 243     bond_inst.each_with_index { |inst, ii|
 244       next if inst.length == 0
 245           btype = ii / 7 + 1
 246           if icode / 10 == 81   #  811 instructions
 247             fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
 248           end
 249       inst.each_with_index { |b, i|
 250         if i % 6 == 0
 251           fp.printf "  %d   %3s", (i >= inst.length - 6 || i % 36 == 30 ? 2 : 1), (i % 36 == 0 ? icode.to_s : "")
 252         end
 253         idx, scode = atoms[b].to_adc
 254         fp.printf "%9d", idx * 100000 + scode
 255         if i % 6 == 5 || i == inst.length - 1
 256           fp.print "\n"
 257           if i == inst.length - 1 || i % 36 == 35
 258             fp.printf "%21s%3d%12s%6.3f\n", "", btype, "", 0.05
 259           end
 260         end
 261       }
 262     }
 263   }
 264
 265   fp.print "  0  1001     0.02\n"  #  Activate hidden line removal
 266   output_bonds.call(821)
 267
 268   #  Atom types
 269   #  Atom type 0=714 (boundary), 1=712 (+principal), 2=716 (+axes), 3=711 (+shades)
 270   atom_inst = Array.new(28) { IntGroup.new }
 271   aattr = (attr ? attr["atoms"] : nil)
 272   atomlist.each { |i|
 273     atype = acol = nil
 274         if aattr
 275           aattr.reverse_each { |at|
 276             if at[0].include?(i)
 277                   atype = at[1] % 4
 278                   acol = (at[2] != 0 ? at[2] - 1 : at[2]) % 7
 279                   break
 280                 end
 281           }
 282         end
 283         if atype == nil
 284           acol = 0
 285       an1 = atoms[i].atomic_number
 286       if an1 == 1
 287         atype = 0
 288       elsif an1 <= 6
 289         atype = 1
 290       else
 291         atype = 3
 292       end
 293         end
 294     idx, scode = atoms[i].to_adc
 295     atom_inst[atype * 7 + acol].add(idx)
 296   }
 297   (atom_inst.count - 1).downto(0) { |ii|
 298     inst = atom_inst[ii]
 299         next if inst.count == 0
 300         fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
 301     i = 0
 302         atype = [714, 712, 716, 711][ii / 7]
 303     while (r = inst.range_at(i)) != nil
 304       fp.printf "  1   %3d\n", atype
 305       fp.printf "%27s%9d%9d\n", "", r.first, r.last
 306       i += 1
 307     end
 308   }
 309
 310   output_bonds.call(811)
 311
 312   #  Close plot
 313   fp.print "      202\n"
 314   fp.print "       -1\n"
 315
 316 end
 317
 318 def savetep(filename)
 319   if natoms == 0
 320         raise MolbyError, "cannot save ORTEP input; the molecule is empty"
 321   end
 322   fp = open(filename, "wb")
 323   export_ortep(fp)
 324   fp.close
 325   return true
 326 end
 327
 328 end
 329
 330 module Math
 331   def sqrt_safe(arg)
 332     arg <= 0.0 ? 0.0 : sqrt(arg)
 333   end
 334 end
 335
 336 #  Best-fit planes
 337 #  Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
 338 #       T. Ito, Acta Cryst 1981, A37, 621-624
 339
 340 #  An object to describe the best-fit plane
 341 #
 342 #  Contains the plane coefficients and the constant (a, b, c, d for ax + by + cz + d = 0),
 343 #  the error matrix, and the metric tensor.
 344 #
 345 class Molby::Plane
 346   attr_accessor :molecule, :group, :coeff, :const, :error_matrix, :metric_tensor
 347   def initialize(mol, group, coeff, const, err, met)
 348     @molecule = mol
 349     @group = group
 350     @coeff = coeff
 351     @const = const
 352     @error_matrix = err
 353     @metric_tensor = met
 354     self
 355   end
 356   def sigma
 357     [sqrt_safe(@error_matrix[0, 0]), sqrt_safe(@error_matrix[1, 1]), sqrt_safe(@error_matrix[2, 2]), sqrt_safe(@error_matrix[3, 3])]
 358   end
 359   def inspect
 360     s = sprintf("Molby::Plane[\n coeff, const = [[%f, %f, %f], %f],\n", @coeff.x, @coeff.y, @coeff.z, @const)
 361     s += sprintf(" sigma = [[%10.4e, %10.4e, %10.4e], %10.4e],\n", *self.sigma)
 362     (0..3).each { |i|
 363       s += (i == 0 ? " error_matrix = [" : "     ")
 364       (0..i).each { |j|
 365         s += sprintf("%12.6e%s", @error_matrix[j, i], (j == i ? (i == 3 ? "],\n" : ",\n") : ","))
 366       }
 367     }
 368     s += sprintf(" molecule = %s\n", @molecule.inspect)
 369     s += sprintf(" group = %s\n", @group.inspect)
 370     (0..3).each { |i|
 371       s += (i == 0 ? " metric_tensor = [" : "     ")
 372       (0..3).each { |j|
 373         s += sprintf("%12.6e%s", @metric_tensor[j, i], (j == 3 ? (i == 3 ? "]]\n" : ",\n") : ","))
 374       }
 375     }
 376     s
 377   end
 378   def distance(ap)
 379     if ap.is_a?(AtomRef)
 380       fr = ap.fract_r
 381       sig = ap.sigma
 382     else
 383       fr = Vector3D[*ap]
 384       sig = Vector3D[0, 0, 0]
 385     end
 386     d = fr.dot(@coeff) + @const
 387     sig1 = (@coeff.x * sig.x) ** 2 + (@coeff.y * sig.y) ** 2 + (@coeff.z * sig.z) ** 2
 388     sig2 = LAMatrix.multiply("t", fr, @error_matrix, fr)[0, 0]
 389     if ap.is_a?(AtomRef) && ap.molecule == @molecule && @group.include?(ap.index)
 390       #  The atom defines the plane
 391       sig0 = sig1 - sig2
 392       sig0 = 0.0 if sig0 < 0.0
 393     else
 394       sig0 = sig1 + sig2
 395     end
 396     return d, sqrt_safe(sig0)
 397   end
 398   def dihedral(plane)
 399     e1 = @error_matrix.submatrix(0, 0, 3, 3)
 400     e2 = plane.error_matrix.submatrix(0, 0, 3, 3)
 401     m = @metric_tensor.submatrix(0, 0, 3, 3)
 402     c = plane.coeff
 403     cos_t = plane.coeff.dot(m * @coeff)
 404     if cos_t > 1.0
 405       cos_t = 1.0
 406     elsif cos_t < -1.0
 407       cos_t = -1.0
 408     end
 409     t = acos(cos_t)
 410     sig_t = (m * e1).trace + (m * e2).trace
 411     if sig_t < t * t
 412       w = 1.0 / sin(t)
 413       sig_t = w * w * (c.dot(LAMatrix.multiply(m, e1, m, c)) + @coeff.dot(LAMatrix.multiply(m, e2, m, @coeff)))
 414     end
 415     t *= 180.0 / PI
 416     sig_t = sqrt_safe(sig_t) * 180.0 / PI
 417     return t, sig_t
 418   end
 419 end
 420
 421 class Molecule
 422
 423 #  Calculate best-fit plane for the given atoms
 424 #  Return value: a Molby::Plane object
 425
 426 def plane(group)
 427
 428   #  Number of atoms
 429   dim = group.length
 430
 431   #  Positional parameters and standard deviations
 432   x = []
 433   sig = []
 434   sig_min = 1e10
 435   each_atom(group) { |ap|
 436     x.push(ap.fract_r)
 437     sig.push(ap.sigma)
 438     if (s = ap.sigma_x) > 0.0 && s < sig_min
 439       sig_min = s
 440     end
 441     if (s = ap.sigma_y) > 0.0 && s < sig_min
 442       sig_min = s
 443     end
 444     if (s = ap.sigma_z) > 0.0 && s < sig_min
 445       sig_min = s
 446     end
 447   }
 448   if sig_min == 1e10
 449     sig_min = 1e-12
 450   end
 451   sig.each { |s|
 452     s.x = sig_min if s.x < sig_min
 453     s.y = sig_min if s.y < sig_min
 454     s.z = sig_min if s.z < sig_min
 455   }
 456
 457   #  The metric tensor of the reciprocal lattice
 458   #  g[j, i] = (ai*).dot(aj*), where ai* and aj* are the reciprocal axis vectors
 459   t = self.cell_transform
 460   if t.nil?
 461     t = Transform.identity
 462   end
 463   g2inv = LAMatrix[t]
 464   g2 = g2inv.inverse
 465   g2[3, 3] = 0.0
 466   g2inv[3, 3] = 0.0
 467   g = LAMatrix.multiply("t", g2, g2)
 468
 469   #  The variance-covariance matrices of the atomic parameters
 470   #  mm[k][n] is a 3x3 matrix describing the correlation between the atoms k and n,
 471   #  and its components are defined as: sigma_k[i] * sigma_n[j] * corr[k, i, n, j],
 472   #  where corr(k, i, n, j) is the correlation coefficients between the atomic parameters
 473   #  k[i] and n[j].
 474   mm = Array.new(dim) { Array.new(dim) }
 475   zero = LAMatrix.zero(3, 3)
 476   dim.times { |k|
 477     dim.times { |n|
 478       mkn = LAMatrix.new(3, 3)
 479       if k == n
 480         3.times { |i|
 481           3.times { |j|
 482             if i == j
 483               mkn[j, i] = sig[k][i] * sig[n][j]
 484             else
 485               #  Inter-coordinate correlation should be implemented here
 486             end
 487           }
 488         }
 489       else
 490         #  Inter-atomic correlation should be implemented here
 491       end
 492       mm[k][n] = (mkn == zero ? zero : mkn)
 493     }
 494   }
 495
 496   #  The variance-covariance matrix of the atom-plance distances
 497   #  m[j, i] = v.transpose * mm[i][j] * v, where v is the plane coefficient vector
 498   #  The inverse of m is the weight matrix
 499   m = LAMatrix.new(dim, dim)
 500
 501   #  The matrix representation of the atomic coordinates
 502   #  y[j, i] = x[i][j] (for j = 0..2), -1 (for j = 3)
 503   #  y * LAMatrix[a, b, c, d] gives the atom-plane distances for each atom
 504   y = LAMatrix.new(4, dim)
 505   dim.times { |i|
 506     y[0, i] = x[i].x
 507     y[1, i] = x[i].y
 508     y[2, i] = x[i].z
 509     y[3, i] = 1.0
 510   }
 511
 512   #  The coefficients to be determined
 513   n0 = LAMatrix[1, 1, 1, 0]
 514   v = LAMatrix[1, 1, 1]     #  The coefficient part
 515
 516   iter = 0
 517   while iter < 20
 518
 519     iter += 1
 520
 521     #  Set zero to the "constant" part, and normalize the "coefficient" part
 522     n0[0, 3] = 0.0
 523     n0 = g2 * n0
 524     n0.multiply!(1.0 / n0.fnorm)
 525     n0 = g2inv * n0
 526     3.times { |i| v[0, i] = n0[0, i] }
 527
 528     #  Build the variance-covariance matrix
 529     dim.times { |i|
 530       dim.times { |j|
 531         m[j, i] = LAMatrix.multiply("t", v, mm[i][j], v)[0, 0]
 532       }
 533     }
 534     c = LAMatrix.multiply("t", y, "i", m, y)
 535
 536     #  Invert c: only the inverse is used in the following, so c is inversed destructively
 537     cinv = c.inverse!
 538
 539     if iter == 1
 540
 541       #  Determine the tentative solution, which is given by the eigenvector of cinv * g
 542       #  for the largest eigenvalue
 543       evals, evecs = (cinv * g).eigenvalues
 544       4.times { |i| n0[0, i] = evecs[3, i] }
 545
 546     else
 547
 548       #  Convert the coefficient vector to the reciprocal space
 549       h = g * n0
 550
 551       #  Determine multiplier
 552       #  In this implementation, the sign of delta-n is opposite from that used in
 553       #  the reference
 554       lam = 1.0 / (LAMatrix.multiply("t", h, cinv, h)[0, 0])
 555
 556       #  Solve the linearized equation
 557       #  (Is the equation 21 in the reference really correct? Shouldn't it read
 558       #   B = 1 - lambda * C.inverse * H* ? )
 559       b = LAMatrix.multiply(lam, cinv, g)
 560       b.sub!(LAMatrix.identity(4))
 561
 562       dn = b * n0
 563       n0 += dn
 564
 565       break if dn[0, 0] ** 2 + dn[0, 1] ** 2 + dn[0, 2] ** 2 < 1e-9
 566
 567     end
 568   end
 569
 570   #  Error matrix = b * cinv * b.transpose
 571   em = LAMatrix.multiply(b, cinv, "t", b)
 572   coeff = Vector3D[n0[0, 0], n0[0, 1], n0[0, 2]]
 573   const = n0[0, 3]
 574
 575   return Molby::Plane.new(self, group, coeff, const, em, g)
 576
 577 end
 578
 579 def cmd_plane
 580   plane_settings = @plane_settings || Hash.new
 581   mol = self
 582   h = Dialog.new("Best-Fit Planes: " + mol.name, "Close", nil) {
 583     refresh_proc = lambda { |it|
 584       n = it[:tag][/\d/].to_i
 585       g = plane_settings["group#{n}"]
 586       if g
 587         str = g.inspect.sub!("IntGroup[", "").sub!("]", "")
 588         set_value("group#{n}", str)
 589         if n == 1 || n == 2
 590           p = mol.plane(g) rescue p = nil
 591           plane_settings["plane#{n}"] = p
 592           if p
 593             coeff = p.coeff
 594             const = p.const
 595             sig = p.sigma
 596             aps = (n == 1 ? "" : "'")
 597             str = sprintf("a%s = %f(%f)\nb%s = %f(%f)\nc%s = %f(%f)\nd%s = %f(%f)",
 598                           aps, coeff.x, sig[0],
 599                           aps, coeff.y, sig[1],
 600                           aps, coeff.z, sig[2],
 601                           aps, const, sig[3])
 602             set_value("result#{n}", str)
 603           else
 604             set_value("result#{n}", "")
 605           end
 606           p1 = plane_settings["plane1"]
 607           p2 = plane_settings["plane2"]
 608           if p1 && p2
 609             t, sig = p1.dihedral(p2)
 610             str = sprintf("%f(%f)", t, sig)
 611             set_value("dihedral", str)
 612           else
 613             set_value("dihedral", "")
 614           end
 615         else
 616           p = plane_settings["plane1"]
 617           if p
 618             str = ""
 619             mol.each_atom(g) { |ap|
 620               d, sig = p.distance(ap)
 621               str += sprintf("%d %f(%f)\n", ap.index, d, sig)
 622             }
 623             str.chomp!
 624           else
 625             str = ""
 626           end
 627           set_value("result#{n}", str)
 628         end
 629       else
 630         set_value("group#{n}", "")
 631         set_value("result#{n}", "")
 632       end
 633     }
 634     set_proc = lambda { |it|
 635       n = it[:tag][/\d/].to_i
 636       sel = mol.selection
 637       if sel.count > 0
 638         str = sel.inspect.sub!("IntGroup[", "").sub!("]", "")
 639         set_value("group#{n}", str)
 640         plane_settings["group#{n}"] = sel
 641       else
 642         plane_settings["group#{n}"] = nil
 643       end
 644       refresh_proc.call(it)
 645     }
 646     text_proc = lambda { |it|
 647       n = it[:tag][/\d/].to_i
 648       str = it[:value].gsub(/[^-.,0-9]/, "")  #  Remove unsane characters
 649       g = eval("IntGroup[#{str}]") rescue g = nil
 650       plane_settings["group#{n}"] = g
 651       refresh_proc.call(it)
 652     }
 653     layout(3,
 654       item(:text, :title=>"Plane 1 (ax + by + cz + d = 0)"),
 655       -1, -1,
 656       item(:text, :title=>"Atoms"),
 657       item(:textfield, :width=>240, :height=>32, :tag=>"group1", :action=>text_proc),
 658       item(:button, :title=>"Set Current Selection", :tag=>"button1", :action=>set_proc),
 659       item(:text, :title=>"Results"),
 660       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result1"),
 661       item(:button, :title=>"Recalculate", :tag=>"refresh1", :action=>refresh_proc),
 662       item(:line),
 663       -1, -1,
 664       item(:text, :title=>"Plane 2 (a'x + b'y + c'z + d' = 0)"),
 665       -1, -1,
 666       item(:text, :title=>"Atoms"),
 667       item(:textfield, :width=>240, :height=>32, :tag=>"group2", :action=>text_proc),
 668       item(:button, :title=>"Set Current Selection", :tag=>"button2", :action=>set_proc),
 669       item(:text, :title=>"Results"),
 670       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result2"),
 671       item(:button, :title=>"Recalculate", :tag=>"refresh2", :action=>refresh_proc),
 672       item(:text, :title=>"Dihedral angle with Plane 1"), -1, -1,
 673       -1,
 674       item(:textfield, :width=>240, :height=>16, :tag=>"dihedral"), -1,
 675       item(:line),
 676       -1, -1,
 677       item(:text, :title=>"Distance from Plane 1"), -1, -1,
 678       item(:text, :title=>"Atoms"),
 679       item(:textfield, :width=>240, :height=>32, :tag=>"group3", :action=>text_proc),
 680       item(:button, :title=>"Set Current Selection", :tag=>"button3", :action=>set_proc),
 681       item(:text, :title=>"Results"),
 682       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result3"),
 683       item(:button, :title=>"Recalculate", :tag=>"refresh3", :action=>refresh_proc)
 684     )
 685     refresh_proc.call(item_with_tag("refresh1"))
 686     refresh_proc.call(item_with_tag("refresh2"))
 687     refresh_proc.call(item_with_tag("refresh3"))
 688         show
 689   }
 690   @plane_settings = plane_settings
 691 end
 692
 693 #  Calculate bond length and angles with standard deviations
 694 #  args can be a single IntGroup (calculate bonds and angles including those atoms)
 695 #  or arrays of 2 or 3 integers (explicitly specifying bonds and angles)
 696 def bond_angle_with_sigma(*args)
 697
 698   if args.length >= 2 || (args.length == 1 && !args[0].is_a?(IntGroup))
 699     #  Bonds and angles are explicitly specified
 700     bonds = []
 701     angles = []
 702     args.each { |arg|
 703       if arg.length == 2 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer)
 704         bonds.push(arg)
 705       elsif arg.length == 3 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer) && arg[2].is_a?(Integer)
 706         angles.push(arg)
 707       else
 708         raise MolbyError, "Invalid argument #{arg.inspect}"
 709       end
 710     }
 711   else
 712     if args.length == 0
 713           g = nil
 714         else
 715           g = args[0]
 716         end
 717         bonds = self.bonds.select { |b|
 718           (g == nil || (g.include?(b[0]) && g.include?(b[1]))) &&
 719             (atoms[b[0]].symop == nil || atoms[b[1]].symop == nil)
 720         }
 721         angles = self.angles.select { |ang|
 722           (g == nil || (g.include?(ang[0]) && g.include?(ang[1]) && g.include?(ang[2]))) &&
 723             (atoms[ang[0]].symop == nil || atoms[ang[1]].symop == nil || atoms[ang[2]].symop == nil)
 724         }
 725   end
 726
 727   #  A list of interatomic distance, its differential coefficients, and other
 728   #  useful quantities.
 729   #  The hash tag is a list [i, j] (i and j are the atom indices) or [i, j, k]
 730   #  (i, j, k are the atom indices, and r(ijk) is the distance between the atom i
 731   #  and the center point between the atoms j and k).
 732   #  The value is a list of following quantities:
 733   #    index 0: rij or r(ijk)
 734   #    index 1-9: d(rij)/d(xij), d(rij)/d(yij), d(rij)/d(zij),
 735   #      d(rij)/da, d(rij)/db, d(rij)/dc, d(rij)/d(alpha), d(rij)/d(beta), d(rij)/d(gamma)
 736   #    index 10: the list of the "base atom"
 737   #    index 11: the list of the transform matrices
 738   dlist = Hash.new
 739
 740   #  A list of fractional coordinates and sigmas
 741   fract = []
 742   sigma = []
 743   each_atom { |ap|
 744     fract.push(ap.fract_r)
 745     sigma.push(ap.sigma)
 746   }
 747
 748   #  A list of base atoms (for symmetry-related atoms) and transform matrices
 749   bases = []
 750   trans = []
 751   symcode = []
 752   trans_i = Transform.identity
 753   each_atom { |ap|
 754     sym = ap.symop
 755     bases.push(sym ? sym[4] : ap.index)
 756         if sym
 757           tr = transform_for_symop(sym).transpose
 758       tr[3, 0] = tr[3, 1] = tr[3, 2] = 0.0
 759         else
 760           tr = trans_i
 761         end
 762     trans.push(tr)
 763         symcode.push(sym ? sprintf("%d_%d%d%d", sym[0] + 1, sym[1] + 5, sym[2] + 5, sym[3] + 5) : ".")
 764   }
 765
 766   #  Unit cell parameter
 767   cell = self.cell
 768   if cell.length == 6
 769     cell.push(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)  #  No sigma information
 770   end
 771   # $a, $b, $c, $alpha, $beta, $gamma, $sig_a, $sig_b, $sig_c, $sig_alpha, $sig_beta, $sig_gamma = self.cell
 772   cos_a = cos(cell[3] * PI / 180.0)
 773   cos_b = cos(cell[4] * PI / 180.0)
 774   cos_c = cos(cell[5] * PI / 180.0)
 775   abc = cell[0] * cell[1] * cos_c
 776   bca = cell[1] * cell[2] * cos_a
 777   cab = cell[2] * cell[0] * cos_b
 778   aa = cell[0] * cell[0]
 779   bb = cell[1] * cell[1]
 780   cc = cell[2] * cell[2]
 781
 782   get_dlist = lambda { |_i, _j, _k|
 783     if _k != nil
 784       if _j > _k
 785         _j, _k = _k, _j
 786       end
 787       _p = dlist[[_i, _j, _k]]
 788       return _p if _p != nil
 789       _p = (dlist[[_i, _j, _k]] = [])
 790       _vij = fract[_i] - (fract[_j] + fract[_k]) * 0.5
 791       _p[10] = [bases[_i], bases[_j], bases[_k]]
 792       _p[11] = [trans[_i], trans[_j], trans[_k]]
 793     else
 794       if _i > _j
 795         _i, _j = _j, _i
 796       end
 797       _p = dlist[[_i, _j]]
 798           return _p if _p != nil
 799       _p = (dlist[[_i, _j]] = [])
 800       _vij = fract[_i] - fract[_j]
 801       _p[10] = [bases[_i], bases[_j]]
 802       _p[11] = [trans[_i], trans[_j]]
 803     end
 804     _xij = _vij.x
 805     _yij = _vij.y
 806     _zij = _vij.z
 807     _dij = sqrt_safe(aa * _xij * _xij + bb * _yij * _yij + cc * _zij * _zij + 2 * bca * _yij * _zij + 2 * cab * _zij * _xij + 2 * abc * _xij * _yij)
 808     _p[0] = _dij
 809     _p[1] = (aa * _xij + abc * _yij + cab * _zij) / _dij
 810     _p[2] = (bb * _yij + bca * _zij + abc * _xij) / _dij
 811     _p[3] = (cc * _zij + cab * _xij + bca * _yij) / _dij
 812     _p[4] = (cell[0] * _xij * _xij + cell[2] * cos_b * _zij * _xij + cell[1] * cos_c * _xij * _yij) / _dij
 813     _p[5] = (cell[1] * _yij * _yij + cell[0] * cos_c * _xij * _yij + cell[2] * cos_a * _yij * _zij) / _dij
 814     _p[6] = (cell[2] * _zij * _zij + cell[1] * cos_a * _yij * _zij + cell[0] * cos_b * _zij * _xij) / _dij
 815     _p[7] = (-cell[1] * cell[2] * sin(cell[3] * PI / 180.0) * _yij * _zij) * (PI / 180.0) / _dij
 816     _p[8] = (-cell[2] * cell[0] * sin(cell[4] * PI / 180.0) * _zij * _xij) * (PI / 180.0) / _dij
 817     _p[9] = (-cell[0] * cell[1] * sin(cell[5] * PI / 180.0) * _xij * _yij) * (PI / 180.0) / _dij
 818     return _p
 819   }
 820
 821   diff_by_rn = lambda { |_dl, _n, _ijk|
 822     #  dl = dlist(i, j)
 823     #  return value: Vector3D[ d(rij)/d(xn), d(rij)/d(yn), d(rij)/d(zn) ]
 824     #  If ijk is true, then dl is dlist(i, j, k)
 825     _dv = Vector3D[_dl[1], _dl[2], _dl[3]]
 826     _c = Vector3D[0, 0, 0]
 827     if _dl[10][0] == _n
 828       _c += _dl[11][0] * _dv
 829     end
 830     if _ijk
 831       if _dl[10][1] == _n
 832         _c -= _dl[11][1] * _dv * 0.5
 833       end
 834       if _dl[10][2] == _n
 835         _c -= _dl[11][2] * _dv * 0.5
 836       end
 837     else
 838       if _dl[10][1] == _n
 839         _c -= _dl[11][1] * _dv
 840       end
 841     end
 842         return _c
 843   }
 844
 845   notate_with_sigma = lambda { |_val, _sig|
 846     if _sig == 0.0
 847       return sprintf "%.4f", _val
 848     end
 849     _lg = log(_sig.abs / 1.95) / log(10.0)
 850     _n = -(_lg.floor)
 851     _val2 = (_val * (10 ** _n) + 0.5).floor * (0.1 ** _n)
 852     return sprintf "%.#{_n}f(%d)", _val2, (_sig * (10 ** _n) + 0.5).floor
 853   }
 854
 855   results = []
 856
 857   c = Vector3D[0, 0, 0]
 858   bonds.each { |b|
 859     i = b[0]
 860     j = b[1]
 861     if i > j
 862       i, j = j, i
 863     end
 864     p = get_dlist.call(i, j, nil)
 865     sig = 0.0
 866     p[10].uniq.each { |k|
 867       s = sigma[k]
 868       next unless s
 869       c = diff_by_rn.call(p, k, false)
 870       #  Test
 871       if nil
 872         apk = atoms[k]
 873         r = apk.fract_r
 874         d0 = calc_bond(i, j)
 875         apk.fract_r = r + Vector3D[0.00001, 0, 0]
 876         amend_by_symmetry(p[10])
 877         d1 = calc_bond(i, j)
 878         apk.fract_r = r + Vector3D[0, 0.00001, 0]
 879         amend_by_symmetry(p[10])
 880         d2 = calc_bond(i, j)
 881         apk.fract_r = r + Vector3D[0, 0, 0.00001]
 882         amend_by_symmetry(p[10])
 883         d3 = calc_bond(i, j)
 884         apk.fract_r = r
 885         amend_by_symmetry(p[10])
 886         printf " dr/dv[%d] cal %10.5g %10.5g %10.5g\n", k, c[0], c[1], c[2]
 887         printf " dr/dv[%d] est %10.5g %10.5g %10.5g\n", k, (d1 - d0) / 0.00001, (d2 - d0) / 0.00001, (d3 - d0) / 0.00001
 888       end
 889       sig += c[0] * c[0] * s[0] * s[0] + c[1] * c[1] * s[1] * s[1] + c[2] * c[2] * s[2] * s[2]
 890     }
 891     6.times { |n|
 892       sig += (p[4 + n] * cell[6 + n]) ** 2
 893     }
 894     sig = sqrt_safe(sig)
 895     results.push([[i, j], notate_with_sigma.call(p[0], sig), symcode[i], symcode[j], nil])
 896   }
 897
 898   angles.each { |ang|
 899     i = ang[0]
 900     j = ang[1]
 901     k = ang[2]
 902     p0 = get_dlist.call(i, j, nil)
 903     p1 = get_dlist.call(j, k, nil)
 904     p2 = get_dlist.call(i, k, nil)
 905     p3 = get_dlist.call(j, i, k)
 906     t123 = acos_safe((p0[0] ** 2 + p1[0] ** 2 - p2[0] ** 2) / (2 * p0[0] * p1[0]))
 907     t124 = acos_safe((p0[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p0[0] * p3[0]))
 908     t324 = acos_safe((p1[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p1[0] * p3[0]))
 909     dtdr12 = -(p0[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p0[0] ** 2) * sin(t124))
 910     dtdr23 = -(p1[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p1[0] ** 2) * sin(t324))
 911     dtdr13 = p2[0] / (sin(t124) * 4 * p0[0] * p3[0]) + p2[0] / (sin(t324) * 4 * p1[0] * p3[0])
 912     dtdr24 = -(p3[0] ** 2 - p0[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * (p3[0] ** 2) * p0[0] * sin(t124)) - (p3[0] ** 2 - p1[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * (p3[0] ** 2) * p1[0] * sin(t324))
 913     pp = (p0[10] + p1[10] + p2[10] + p3[10]).uniq
 914     sig = 0.0
 915     pp.each { |n|
 916       s = sigma[n]
 917       next unless s
 918       c = Vector3D[0, 0, 0]
 919       c += diff_by_rn.call(p0, n, false) * (dtdr12 * 180.0 / PI)
 920       c += diff_by_rn.call(p1, n, false) * (dtdr23 * 180.0 / PI)
 921       c += diff_by_rn.call(p2, n, false) * (dtdr13 * 180.0 / PI)
 922       c += diff_by_rn.call(p3, n, true) * (dtdr24 * 180.0 / PI)
 923       sig += c[0] * c[0] * s[0] * s[0] + c[1] * c[1] * s[1] * s[1] + c[2] * c[2] * s[2] * s[2]
 924     }
 925     dd = dtdr12 * p0[4] + dtdr23 * p1[4] + dtdr13 * p2[4] + dtdr24 * p3[4]
 926     sig += dd * dd * cell[6] * cell[6]
 927     dd = dtdr12 * p0[5] + dtdr23 * p1[5] + dtdr13 * p2[5] + dtdr24 * p3[5]
 928     sig += dd * dd * cell[7] * cell[7]
 929     dd = dtdr12 * p0[6] + dtdr23 * p1[6] + dtdr13 * p2[6] + dtdr24 * p3[6]
 930     sig += dd * dd * cell[8] * cell[8]
 931     dd = dtdr12 * p0[7] + dtdr23 * p1[7] + dtdr13 * p2[7] + dtdr24 * p3[7]
 932     sig += dd * dd * cell[9] * cell[9]
 933     dd = dtdr12 * p0[8] + dtdr23 * p1[8] + dtdr13 * p2[8] + dtdr24 * p3[8]
 934     sig += dd * dd * cell[10] * cell[10]
 935     dd = dtdr12 * p0[9] + dtdr23 * p1[9] + dtdr13 * p2[9] + dtdr24 * p3[9]
 936     sig += dd * dd * cell[11] * cell[11]
 937     sig = sqrt_safe(sig)
 938     results.push([[i, j, k], notate_with_sigma.call(t123*180/PI, sig), symcode[i], symcode[j], symcode[k]])
 939   }
 940   results
 941 end
 942
 943 def cmd_bond_angle_with_sigma
 944   mol = self
 945   Dialog.new("Bond & Angle with Sigma:" + self.name, nil, nil) {
 946     values = []
 947         clicked = []
 948         sel = mol.selection
 949         on_get_value = lambda { |it, row, col| values[row][col + 2] }
 950         atom_name = lambda { |ap|
 951           (ap.molecule.nresidues >= 2 ? "#{ap.res_seq}:" : "") + ap.name
 952         }
 953     layout(1,
 954           item(:table, :width=>480, :height=>300, :tag=>"table",
 955             :columns=>[["Bond/Angle", 140], "value(sigma)", "symop1", "symop2", "symop3"],
 956                 :on_count=> lambda { |it| values.count },
 957                 :on_get_value=> on_get_value),
 958           layout(2,
 959             item(:view, :width=>480, :height=>1), -1,
 960             item(:button, :title=>"Export to Clipboard", :tag=>"dump",
 961                   :action=>lambda { |item|
 962                     s = ""; values.each { |val| s += val[2..-1].join("  ") + "\n" }; export_to_clipboard(s)
 963                   },
 964                   :enabled=>false),
 965                 [item(:button, :title=>"Close", :action=>lambda { |item| hide }), {:align=>:right}])
 966         )
 967         on_document_modified = lambda { |*d|
 968           puts "on_document_modified called: newsel = #{mol.selection}, sel = #{sel}"
 969           newsel = mol.selection
 970           val1 = val2 = nil
 971           if sel != newsel
 972             n1 = newsel.count
 973                 if n1 == 0
 974                   #  Clear selection
 975                   clicked.clear
 976                 elsif n1 == 1 && !clicked.include?(newsel[0])
 977                   #  New atom
 978                   clicked.clear
 979                   clicked[0] = newsel[0]
 980                   val1 = [[clicked[0]], "---", "", "", ""]
 981                 elsif n1 == 2 && clicked.length == 1
 982                   #  New bond
 983                   if newsel[0] == clicked[0]
 984                     clicked[1] = newsel[1]
 985                   elsif newsel[1] == clicked[0]
 986                     clicked[1] = newsel[0]
 987                   else
 988                     clicked[0] = newsel[0]
 989                         clicked[1] = newsel[1]
 990                   end
 991                   val1 = mol.bond_angle_with_sigma(clicked)[0]
 992                 elsif n1 == 3 && clicked.length == 2 && newsel.include?(clicked[0]) && newsel.include?(clicked[1])
 993                   #  New angle
 994                   clicked[2] = (newsel - clicked)[0]
 995                   val1 = mol.bond_angle_with_sigma(clicked[1..2])[0]
 996                   val2 = mol.bond_angle_with_sigma(clicked)[0]
 997                 else
 998                   return
 999                 end
1000                 [val1, val2].each { |val|
1001                   next unless val
1002                   label = ""
1003                   n = val[0].length
1004                   n.times { |i|
1005                     label += atom_name.call(mol.atoms[val[0][i]])
1006                         if val[0][i + 1]
1007                           label += (mol.bond_exist?(val[0][i], val[0][i + 1]) ? "-" : "...")
1008                         end
1009                   }
1010                   val[1, 0] = label
1011                   val.unshift(n)  #  val = [count, indices, label, value(sigma), symop1, symop2, symop3]
1012                 }
1013             if values[-1] && values[-1][0] == 1
1014                   values.pop
1015                 end
1016                 values.push(val1) if val1
1017                 values.push(val2) if val2
1018                 sel = newsel
1019                 item_with_tag("table")[:refresh] = true
1020                 item_with_tag("dump")[:enabled] = (values.length > 0)
1021           end
1022         }
1023     listen(mol, "documentModified", on_document_modified)
1024         listen(mol, "documentWillClose", lambda { |*d| hide } )
1025         on_document_modified.call
1026         show
1027   }
1028 end
1029
1030 def find_expanded_atom(base, symop)
1031   return base unless symop
1032   symopb = symop + [base]
1033   self.each_atom { |ap|
1034     if ap.symop == symopb
1035       return ap.index
1036     end
1037   }
1038   nil
1039 end
1040
1041 def complete_by_symmetry
1042   if self.box == nil
1043     raise "Unit cell should be given"
1044   end
1045   verbose = false
1046   avec, bvec, cvec = self.box
1047   syms = []
1048   self.nsymmetries.times { |n|
1049     syms.push(transform_for_symop([n, 0, 0, 0], true))
1050   }
1051   frags = []
1052   self.each_fragment { |f|
1053     frags.push(f)
1054   }
1055   close_pairs = []
1056   self.each_atom { |ap|
1057     #  Find if ap is contained in expansion
1058         next if ap.symop != nil  #  Already expanded
1059     rad = Parameter.builtin.elements[ap.atomic_number].radius
1060     syms.each_with_index { |sym, sym_idx|
1061       27.times { |n|
1062         dx = n % 3 - 1
1063         dy = (n / 3) % 3 - 1
1064         dz = (n / 9) % 3 - 1
1065         next if dx == 0 && dy == 0 && dz == 0 && sym_idx == 0
1066         symop = [sym_idx, dx, dy, dz]
1067                 next if self.find_expanded_atom(ap.index, symop)
1068         r = sym * ap.r + avec * dx + bvec * dy + cvec * dz
1069         close_atoms = self.find_close_atoms(r, 1.2, rad)
1070         if close_atoms.length > 0
1071           #  ap * [sym, dx, dy, dz] is included in the expansion
1072           close_atoms.each { |ca|
1073             next if ca > ap.index
1074             i1 = frags.index { |f| f.include?(ca) }
1075             i2 = frags.index { |f| f.include?(ap.index) }
1076             pp = close_pairs.find { |p1| p1[0] == i1 && p1[1] == i2 && p1[2] == symop }
1077                         if pp == nil
1078                           pp = [i1, i2, symop, [], [], [], []]
1079                           close_pairs.push(pp)
1080                     end
1081             if (r - atoms[ca].r).length2 > 0.0001
1082                           #  Normal case (bond between ca and ap)
1083                           pp[3].push(ca)
1084                           pp[4].push(ap.index)
1085                         else
1086               #  Special position
1087               pp[5].push(ca)
1088               pp[6].push(ap.index)
1089             end
1090           }
1091         end
1092       }
1093     }
1094   }
1095   puts "close_pairs = #{close_pairs}" if verbose
1096   expand_pairs = lambda { |i1, symop, depth|
1097     #  Find expanded fragment that is close to fragment [i1, symop]
1098     next [] if depth > 16
1099     retval = []
1100     close_pairs.each { |pp|
1101       next if i1 != nil && pp[0] != i1
1102           next if pp[3].count == 0   #  No real expansion
1103       #  Multiply two symops
1104       if symop
1105         symop2 = symop_for_transform(transform_for_symop(pp[2]) * transform_for_symop(symop))
1106         puts "multiply two symops: #{pp[2].inspect} * #{symop.inspect} -> #{symop2.inspect}" if verbose
1107       else
1108         symop2 = pp[2]
1109       end
1110       if symop2[0] == 0
1111         #  Translation only
1112         if symop2[1] == 0 && symop2[2] == 0 && symop2[3] == 0
1113                   #  No expansion, but pp[3] and pp[4] are still needed to make bonds
1114                   symop2 = nil
1115             end
1116                 #  Expand the atoms only at the marginal
1117             retval.push([nil, symop, symop2, pp[3], pp[4], pp[5], pp[6]])
1118       else
1119         #  Expand fragment
1120         retval.push([pp[1], symop, symop2, pp[3], pp[4], pp[5], pp[6]])
1121         retval.concat(expand_pairs.call(pp[1], symop2, depth + 1))
1122       end
1123     }
1124     next retval
1125   }
1126   ex_pairs = expand_pairs.call(nil, nil, 0)
1127   puts "ex_pairs = #{ex_pairs}" if verbose
1128
1129   #  Expand fragments
1130   duplicates = []
1131   ex_pairs.each { |ex|
1132     #  ex[0]: fragment to expand, ex[1], ex[2]: symop,
1133     #  ex[3], ex[4]: atoms to make bonds, ex[5], ex[6]: duplicate atoms
1134     n1 = self.natoms
1135     if ex[0] == nil
1136       f = IntGroup[ex[4]]
1137     else
1138       f = frags[ex[0]]
1139     end
1140         if ex[2] != nil
1141       self.expand_by_symmetry(f, ex[2][0], ex[2][1], ex[2][2], ex[2][3], true)
1142       puts "expand(#{f}, #{ex[2]}) -> atoms[#{n1}..#{self.natoms - n1}]" if verbose
1143         end
1144     ex[3].each_with_index { |x, i|
1145       x1 = self.find_expanded_atom(x, ex[1])
1146           if ex[2] != nil
1147         x2 = f.index(ex[4][i]) + n1  #  New index of the expanded atom
1148           else
1149             x2 = ex[4][i]   #  Base atom
1150           end
1151       create_bond(x1, x2)
1152       puts "create_bond(#{x1}, #{x2})" if verbose
1153     }
1154     #  Register duplicate atoms
1155         if ex[2] != nil
1156       ex[5].each_with_index { |x, i|
1157         x1 = self.find_expanded_atom(x, ex[1])
1158         x2 = f.index(ex[6][i]) + n1
1159         duplicates.each { |d|
1160           if d.include?(x1)
1161             d.push(x2)
1162             x2 = nil
1163             break
1164           end
1165         }
1166         if x2
1167           duplicates.push([x1, x2])
1168         end
1169       }
1170     end
1171   }
1172   puts "duplicates = #{duplicates}" if verbose
1173
1174   #  Remove duplicate atoms
1175   rem = []
1176   duplicates.each { |d|
1177     d.each_with_index { |dn, i|
1178       next if i == 0
1179       #  Remake bonds
1180       self.atoms[dn].connects.each { |nn|
1181         create_bond(d[0], nn)
1182         puts "create_bond(#{d[0]}, #{nn})" if verbose
1183       }
1184       rem.push(dn)
1185     }
1186   }
1187   remove(rem)
1188   puts "remove(#{rem})" if verbose
1189
1190 end
1191
1192 def create_packing_diagram
1193   if self.box == nil
1194     error_message_box "Unit cell is not defined."
1195         return
1196   end
1197   expansion_box = (@expansion_box ||= [0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
1198   h = Dialog.run("Create Packing Diagram") {
1199     layout(4,
1200           item(:text, :title=>"Specify the expansion limits for each axis:"),
1201           -1, -1, -1,
1202           item(:text, :title=>"x"),
1203           item(:textfield, :width=>80, :tag=>"xmin", :value=>sprintf("%.1f", expansion_box[0].to_f)),
1204           item(:text, :title=>"..."),
1205           item(:textfield, :width=>80, :tag=>"xmax", :value=>sprintf("%.1f", expansion_box[1].to_f)),
1206           item(:text, :title=>"y"),
1207           item(:textfield, :width=>80, :tag=>"ymin", :value=>sprintf("%.1f", expansion_box[2].to_f)),
1208           item(:text, :title=>"..."),
1209           item(:textfield, :width=>80, :tag=>"ymax", :value=>sprintf("%.1f", expansion_box[3].to_f)),
1210           item(:text, :title=>"z"),
1211           item(:textfield, :width=>80, :tag=>"zmin", :value=>sprintf("%.1f", expansion_box[4].to_f)),
1212           item(:text, :title=>"..."),
1213           item(:textfield, :width=>80, :tag=>"zmax", :value=>sprintf("%.1f", expansion_box[5].to_f)))
1214   }
1215   if h[:status] == 0
1216     @expansion_box[0] = h["xmin"].to_f
1217     @expansion_box[1] = h["xmax"].to_f
1218     @expansion_box[2] = h["ymin"].to_f
1219     @expansion_box[3] = h["ymax"].to_f
1220     @expansion_box[4] = h["zmin"].to_f
1221     @expansion_box[5] = h["zmax"].to_f
1222   else
1223     return
1224   end
1225
1226   #  Enumerate all atoms within the expansion box
1227   syms = self.symmetries
1228   h = Hash.new
1229   xmin = @expansion_box[0]
1230   xmax = @expansion_box[1]
1231   ymin = @expansion_box[2]
1232   ymax = @expansion_box[3]
1233   zmin = @expansion_box[4]
1234   zmax = @expansion_box[5]
1235   xmin_d = xmin.floor
1236   xmax_d = xmax.floor - 1
1237   ymin_d = ymin.floor
1238   ymax_d = ymax.floor - 1
1239   zmin_d = zmin.floor
1240   zmax_d = zmax.floor - 1
1241   syms.each_with_index { |sym, i|
1242     each_atom { |ap|
1243       fr = sym * ap.fract_r
1244           dx = fr.x.floor
1245           dy = fr.y.floor
1246           dz = fr.z.floor
1247           fr.x -= dx
1248           fr.y -= dy
1249           fr.z -= dz
1250           symopi = i * 1000000 + (50 - dx) * 10000 + (50 - dy) * 100 + 50 - dz
1251           (h[symopi] ||= []).push(ap.index)
1252           xmin_d.upto(xmax_d) { |xd|
1253             next if fr.x + xd < xmin || fr.x + xd > xmax
1254             ymin_d.upto(ymax_d) { |yd|
1255               next if fr.y + yd < ymin || fr.y + yd > ymax
1256                   zmin_d.upto(zmax_d) { |zd|
1257                     next if fr.z + zd < zmin || fr.z + zd > zmax
1258                         next if xd == 0 && yd == 0 && zd == 0
1259                     symopid = symopi + xd * 10000 + yd * 100 + zd
1260                         (h[symopid] ||= []).push(ap.index)
1261                   }
1262                 }
1263           }
1264         }
1265   }
1266   h.keys.sort.each { |key|
1267     sym = key / 1000000
1268         dx = key % 1000000 / 10000 - 50
1269         dy = key % 10000 / 100 - 50
1270         dz = key % 100 - 50
1271         next if sym == 0 && dx == 0 && dy == 0 && dz == 0
1272 #       puts "[#{sym},#{dx},#{dy},#{dz}]: #{h[key].inspect}"
1273         expand_by_symmetry(IntGroup[h[key]], sym, dx, dy, dz)
1274   }
1275 end
1276
1277 def cmd_show_ortep
1278   mol = self
1279   tmp = create_temp_dir("ortep", mol.name)
1280   tepexe = "#{ResourcePath}/ortep3/ortep3"
1281   if $platform == "win"
1282     tepexe += ".exe"
1283   end
1284   tepdata = []
1285   tepbounds = [0, 0, 400, 400]
1286   descs = {
1287         "Atoms"=>[
1288           ["all", "Shaded", "black", ""],
1289           ["Li-Ar", "Principals", "black", ""],
1290           ["C", "Boundary", "black", ""],
1291           ["H", "Fixed", "black", "0.1"]
1292         ],
1293         "Bonds"=>[
1294           ["all", "all", "4 Shades", "black"],
1295           ["Li-Ar", "Li-Ar", "3 Shades", "black"],
1296           ["C", "Li-Ar", "2 Shades", "black"],
1297           ["H", "all", "Open", "black"]
1298         ]
1299   }
1300   Dialog.new("Show ORTEP:" + mol.name, nil, nil, :resizable=>true, :has_close_box=>true) {
1301     tepview = nil  #  Forward declaration
1302         tab = "Atoms"
1303         columns = {
1304           "Atoms"=>[["atom list", 65], ["type", 80], ["color", 40], ["radius", 50]],
1305           "Bonds"=>[["list1", 65], ["list2", 65], ["type", 80], ["color", 40]]
1306         }
1307         atom_types = ["Boundary", "Principals", "Axes", "Shaded", "Fixed"]
1308         bond_types = ["Open", "1 Shade", "2 Shades", "3 Shades", "4 Shades"]
1309         colors = ["black", "red", "green", "blue", "cyan", "magenta", "yellow"]
1310         #  ORTEP attributes
1311         get_ortep_attr = lambda {
1312           attr = Hash.new
1313           make_atom_list = lambda { |s, ln|
1314             ss = s.scan(/[-.0-9A-Za-z]+/)
1315                 puts ss.inspect
1316                 ss.inject(IntGroup[]) { |r, it|
1317                   if it == "all"
1318                     r.add(mol.all)
1319                         next r
1320                   elsif it =~ /-|(\.\.)/
1321                     s0 = Regexp.last_match.pre_match
1322                         s1 = Regexp.last_match.post_match
1323                   else
1324                     s0 = s1 = it
1325                   end
1326                   if s0 =~ /^[0-9]+$/
1327                     #  Atomic numbers
1328                         s0 = s0.to_i
1329                         if s1 !~ /^[0-9]+$/ || (s1 = s1.to_i) < s0
1330                           raise "Bad atom list specification: #{s} in line #{ln}"
1331                         end
1332                         r.add(s0..s1)
1333                   else
1334                     #  Atomic symbols
1335                         s0 = Parameter.builtin.elements.find_index { |p| p.name == s0.capitalize }
1336                         s1 = Parameter.builtin.elements.find_index { |p| p.name == s1.capitalize }
1337                         if s0 == nil || s1 == nil
1338                           raise "Bad element symbol specification: #{s} in line #{ln}"
1339                         end
1340                         (s0..s1).each { |an|
1341                           r.add(mol.atom_group { |ap| ap.atomic_number == an } )
1342                         }
1343                   end
1344                   r
1345                 }
1346       }
1347           #  Atoms
1348           at = []
1349           descs["Atoms"].each_with_index { |d, idx|
1350                 list = make_atom_list.call(d[0], idx)
1351                 type = atom_types.find_index { |it| it == d[1] } || 0
1352                 col = (colors.find_index { |it| it == d[2] } || 0) + 1
1353                 rad = d[3]
1354                 at.push([list, type, col, rad])
1355           }
1356           attr["atoms"] = at
1357           #  Bonds
1358           bd = []
1359           descs["Bonds"].each_with_index { |d, idx|
1360                 list1 = make_atom_list.call(d[0], idx)
1361                 list2 = make_atom_list.call(d[1], idx)
1362                 type = bond_types.find_index { |it| it == d[2] } || 0
1363                 col = (colors.find_index { |it| it == d[3] } || 0) + 1
1364                 bd.push([list1, list2, type, col])
1365           }
1366           attr["bonds"] = bd
1367           attr
1368         }
1369     on_update_ortep = lambda { |it|
1370           #  ORTEP attributes
1371           attr = get_ortep_attr.call
1372           #  Create ORTEP input in the temporary directory
1373           open(tmp + "/TEP.IN", "w") { |fp| mol.export_ortep(fp, attr) }
1374           #  Run ORTEP
1375           cwd = Dir.pwd
1376           Dir.chdir(tmp)
1377           if FileTest.exist?("TEP001.PRN")
1378             File.unlink("TEP001.PRN")
1379           end
1380           pid = call_subprocess(tepexe, "Running ORTEP", nil, "NUL", "NUL")
1381           if FileTest.exist?("TEP001.PRN")
1382             File.rename("TEP001.PRN", "TEP001.ps")
1383           end
1384           Dir.chdir(cwd)
1385           if pid != 0
1386             msg = "ORTEP execution in #{tmp} failed with status #{pid}."
1387             message_box(msg, "ORTEP Failed", :ok, :warning)
1388           else
1389             open(tmp + "/TEP001.ps", "r") { |fp|
1390                   tepdata.clear
1391                   tepbounds = [100000, 100000, -100000, -100000]
1392                   fp.each_line { |ln|
1393                     ln.chomp!
1394                     x, y, c = ln.split
1395                         if ln =~ /setrgbcolor/
1396                           tepdata.push(["color", x.to_f, y.to_f, c.to_f])
1397                           next
1398                         elsif ln =~ /setgray/
1399                           x = x.to_f
1400                           tepdata.push(["color", x, x, x])
1401                           next
1402                         end
1403                         x = x.to_i
1404                         y = y.to_i
1405                         if c == "m"
1406                           tepdata.push([x, y])
1407                         elsif c == "l"
1408                           tepdata[-1].push(x, y)
1409                         else
1410                           next
1411                         end
1412                         tepbounds[0] = x if x < tepbounds[0]
1413                         tepbounds[1] = y if y < tepbounds[1]
1414                         tepbounds[2] = x if x > tepbounds[2]
1415                         tepbounds[3] = y if y > tepbounds[3]
1416                   }
1417                   #  [x, y, width, height]
1418                   tepbounds[2] -= tepbounds[0]
1419                   tepbounds[3] -= tepbounds[1]
1420                 }
1421                 tepview.refresh_rect(tepview[:frame])
1422           end
1423         }
1424         on_draw_ortep = lambda { |it|
1425           clear
1426           frame = it[:frame]
1427           brush(:color=>[1, 1, 1, 1])
1428           pen(:color=>[1, 1, 1, 1])
1429           draw_rectangle(frame)
1430           pen(:color=>[0, 0, 0, 1])
1431           rx = (frame[2] - 10.0) / tepbounds[2]
1432           ry = (frame[3] - 10.0) / tepbounds[3]
1433           rx = ry if rx > ry
1434           dx = (frame[2] - tepbounds[2] * rx) * 0.5
1435           dy = (frame[3] + tepbounds[3] * rx) * 0.5
1436           tepdata.each { |d|
1437             if d[0] == "color"
1438                   pen(:color=>[d[1], d[2], d[3], 1])
1439                   next
1440                 end
1441             x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
1442                 y0 = (dy - (d[1] - tepbounds[1]) * rx).to_i
1443                 (d.length / 2 - 1).times { |i|
1444                   x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
1445                   y1 = (dy - (d[i * 2 + 3] - tepbounds[1]) * rx).to_i
1446                   draw_line(x0, y0, x1, y1)
1447                   x0 = x1
1448                   y0 = y1
1449                 }
1450           }
1451         }
1452         on_export_eps = lambda { |fname|
1453           frame = item_with_tag("ortep")[:frame].dup
1454           dx = dy = 5
1455           rx = (frame[2] - dx * 2) / tepbounds[2]
1456           ry = (frame[3] - dy * 2) / tepbounds[3]
1457           rx = ry if rx > ry
1458           frame[2] = (rx * tepbounds[2] + dx * 2).to_i
1459           frame[3] = (rx * tepbounds[3] + dy * 2).to_i
1460           #  Assumes A4 paper and center the bounding box
1461           frame[0] = (595 - frame[2]) / 2
1462           frame[1] = (842 - frame[3]) / 2
1463           xmax = frame[0] + frame[2]
1464           ymax = frame[1] + frame[3]
1465           open(fname, "w") { |fp|
1466             fp.print "%!PS-Adobe-3.0 EPSF-3.0
1467 %%Creator: Molby
1468 %%BoundingBox: #{frame[0]} #{frame[1]} #{xmax} #{ymax}
1469 %%Pages: 1
1470 %%Orientation: Landscape
1471 %%BeginProlog
1472 /m {moveto} def
1473 /l {lineto} def
1474 %%EndProlog
1475 %%Page: 1 1
1476 %%BeginPageSetup
1477 0 setgray 1 setlinecap 1 setlinewidth
1478 %%EndPageSetup
1479 "
1480             tepdata.each { |d|
1481               if d[0] == "color"
1482                     fp.printf "%.3f %.3f %.3f setrgbcolor\n", d[1], d[2], d[3]
1483                     next
1484                   end
1485               x0 = frame[0] + dx + (d[0] - tepbounds[0]) * rx
1486                   y0 = frame[1] + dy + (d[1] - tepbounds[1]) * rx
1487                   fp.printf "%.2f %.2f m\n", x0, y0
1488                   (d.length / 2 - 1).times { |i|
1489                     x1 = frame[0] + dx + (d[i * 2 + 2] - tepbounds[0]) * rx
1490                     y1 = frame[1] + dy + (d[i * 2 + 3] - tepbounds[1]) * rx
1491                     fp.printf "%.2f %.2f l\n", x1, y1
1492                     x0 = x1
1493                     y0 = y1
1494                   }
1495                   fp.print "stroke\n"
1496             }
1497         fp.print "showpage\n%%Trailer\n%%EOF\n"
1498           }
1499         }
1500         on_export_bitmap = lambda { |fname|
1501           frame = item_with_tag("ortep")[:frame].dup
1502           dx = dy = 5
1503           scale = 4.0
1504           rx = (frame[2] - dx * 2) * scale / tepbounds[2]
1505           ry = (frame[3] - dy * 2) * scale / tepbounds[3]
1506           rx = ry if rx > ry
1507           frame[2] = (rx * tepbounds[2] + dx * 2).to_i
1508           frame[3] = (rx * tepbounds[3] + dy * 2).to_i
1509       bmp = Bitmap.new(frame[2], frame[3], 32)
1510           bmp.focus_exec {
1511             clear
1512             pen(:color=>[0, 0, 0, 1], :width=>scale)
1513             tepdata.each { |d|
1514               if d[0] == "color"
1515                     pen(:color=>[d[1], d[2], d[3], 1])
1516                     next
1517                   end
1518               x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
1519                   y0 = (dy + (tepbounds[3] - d[1] + tepbounds[1]) * rx).to_i
1520                   (d.length / 2 - 1).times { |i|
1521                     x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
1522                     y1 = (dy + (tepbounds[3] - d[i * 2 + 3] + tepbounds[1]) * rx).to_i
1523                     draw_line(x0, y0, x1, y1)
1524                     x0 = x1
1525                     y0 = y1
1526                   }
1527             }
1528           }
1529           bmp.save_to_file(fname)
1530         }
1531         on_export = lambda { |it|
1532           basename = (mol.path ? File.basename(mol.path, ".*") : mol.name)
1533       fname = Dialog.save_panel("Export ORTEP file:", mol.dir, basename + ".eps",
1534             "Encapsulated PostScript (*.eps)|*.eps|PNG (*.png)|*.png|TIFF (*.tif)|*.tif|ORTEP Instruction (*.tep)|*.tep|All files|*.*")
1535           return if !fname
1536           ext = File.extname(fname).downcase
1537           if ext == ".eps"
1538             on_export_eps.call(fname)
1539           elsif ext == ".png" || ext == ".tif" || ext == ".tiff"
1540             on_export_bitmap.call(fname)
1541           elsif ext == ".tep"
1542             filecopy("#{tmp}/TEP.IN", fname)
1543       end
1544         }
1545         on_document_modified = lambda { |m|
1546         }
1547         #  Close handler (called when the close box is pressed or the document is closed)
1548         @on_close = lambda { |*d|
1549           cleanup_temp_dir(tmp)
1550           tmp = nil
1551           true
1552         }
1553         on_select_tab = lambda { |it|
1554           tab = ["Atoms", "Bonds"][it[:value]]
1555           puts "tab = #{tab}"
1556           table = item_with_tag("table")
1557           table[:columns] = columns[tab]
1558           table[:selection] = IntGroup[]
1559           table[:refresh] = true
1560         }
1561         get_count = lambda { |it| descs[tab].count }
1562         get_value = lambda { |it, row, col| descs[tab][row][col] }
1563         has_popup_menu = lambda { |it, row, col|
1564           if tab == "Atoms"
1565             if col == 1
1566               return atom_types
1567                 elsif col == 2
1568                   return colors
1569                 end
1570           elsif tab == "Bonds"
1571                 if col == 2
1572               return bond_types
1573                 elsif col == 3
1574                   return colors
1575                 end
1576           end
1577           return nil
1578         }
1579         on_selection_changed = lambda { |it|
1580           sel = it[:selection]
1581           item_with_tag("remove")[:enabled] = (sel == nil || sel.count == 0 ? false : true)
1582         }
1583         is_item_editable = lambda { |it, row, col|
1584           return has_popup_menu.call(it, row, col) == nil
1585         }
1586         popup_menu_selected = lambda { |it, row, col, sel|
1587           titles = has_popup_menu.call(it, row, col)
1588           return nil if titles == nil
1589           descs[tab][row][col] = titles[sel]
1590           it[:refresh] = true
1591         }
1592         set_value = lambda { |it, row, col, val|
1593           descs[tab][row][col] = val
1594           it[:refresh] = true
1595         }
1596         add_to_table = lambda { |it|
1597           table = item_with_tag("table")
1598           d = descs[tab][-1].dup
1599           descs[tab].push(d)
1600           table[:refresh] = true
1601           table[:selection] = IntGroup[descs[tab].count - 1]
1602           puts descs[tab].inspect
1603         }
1604         remove_from_table = lambda { |it|
1605           table = item_with_tag("table")
1606           sel = table[:selection]
1607           sel.reverse_each { |n|
1608             descs[tab][n, 1] = []
1609           }
1610           table[:refresh] = true
1611         }
1612         layout(2,
1613           layout(1,
1614             item(:popup, :subitems=>["Atoms", "Bonds"], :action=>on_select_tab),
1615                 item(:table, :width=>240, :height=>380, :flex=>[0,0,0,0,1,1], :tag=>"table",
1616                   :columns=>columns["Atoms"],
1617                   :on_count=>get_count,
1618                   :on_get_value=>get_value,
1619                   :on_set_value=>set_value,
1620                   :on_selection_changed=>on_selection_changed,
1621                   :is_item_editable=>is_item_editable,
1622                   :has_popup_menu=>has_popup_menu,
1623                   :on_popup_menu_selected=>popup_menu_selected),
1624                 layout(2,
1625                   item(:button, :title=>"Add", :tag=>"add", :action=>add_to_table),
1626                   item(:button, :title=>"Remove", :tag=>"remove", :action=>remove_from_table, :enabled=>false),
1627                   :flex=>[0,1,0,0,0,0]),
1628             :flex=>[0,0,0,0,0,1]),
1629           layout(1,
1630             item(:view, :width=>400, :height=>400, :tag=>"ortep", :on_paint=>on_draw_ortep, :flex=>[0,0,0,0,1,1]),
1631             layout(2,
1632               item(:button, :title=>"Refresh", :action=>on_update_ortep, :align=>:left),
1633                   item(:button, :title=>"Export as...", :action=>on_export, :align=>:right),
1634                   :flex=>[0,1,0,0,0,0]),
1635             :flex=>[0,0,0,0,1,1]),
1636           :flex=>[0,0,0,0,1,1]
1637         )
1638         set_min_size(480, 200)
1639         tepview = item_with_tag("ortep")
1640         listen(mol, "documentModified", on_document_modified)
1641         listen(mol, "documentWillClose", lambda { |m| close } )
1642         on_document_modified.call(nil)
1643     on_update_ortep.call(nil)
1644         show
1645   }
1646 end
1647
1648 if lookup_menu("Best-fit Planes...") < 0
1649   register_menu("", "")
1650   register_menu("Best-fit Planes...", :cmd_plane)
1651   register_menu("Bonds and Angles with Sigma...", :cmd_bond_angle_with_sigma)
1652   register_menu("Complete by Symmetry", :complete_by_symmetry)
1653   register_menu("Create Packing Diagram...", :create_packing_diagram)
1654   register_menu("Show ORTEP", :cmd_show_ortep)
1655 end
1656
1657 end