Scripts/crystal.rb

   1 # coding: utf-8
   2 #
   3 #  crystal.rb
   4 #
   5 #  Created by Toshi Nagata.
   6 #  Copyright 2012 Toshi Nagata. All rights reserved.
   7 #
   8 # This program is free software; you can redistribute it and/or modify
   9 # it under the terms of the GNU General Public License as published by
  10 # the Free Software Foundation version 2 of the License.
  11 #
  12 # This program is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16
  17 #  Definition for use in ORTEP
  18 class AtomRef
  19   def to_adc
  20     sym = self.symop
  21     if sym == nil
  22       idx = self.index + 1
  23       symcode = 55501
  24     else
  25       idx = sym[4] + 1
  26       symcode = (sym[1] + 5) * 10000 + (sym[2] + 5) * 1000 + (sym[3] + 5) * 100 + sym[0] + 1
  27     end
  28     return idx, symcode
  29   end
  30 end
  31
  32 module Kernel
  33
  34 def symmetry_to_string(sym)
  35   #  Sym is a Transform object defined in crystallographic coordinates
  36   str = ""
  37   3.times { |j|
  38     s = ""
  39         a = sym[3, j]
  40         if a.abs > 1e-4
  41           if a > 0.0 && s != ""
  42             s += "+"
  43           elsif a < 0.0
  44             s += "-"
  45           end
  46           a = a.abs
  47           if (a - (n = (a + 0.5).floor)).abs < 1e-4
  48             #  integer
  49                 s += n.to_s
  50           elsif ((a * 2) - (n = (a * 2 + 0.5).floor)).abs < 1e-4
  51             #  n/2
  52             s += n.to_s + "/2"
  53           elsif ((a * 3) - (n = (a * 3 + 0.5).floor)).abs < 1e-4
  54             #  n/3
  55                 s += n.to_s + "/3"
  56           elsif ((a * 4) - (n = (a * 4 + 0.5).floor)).abs < 1e-4
  57             #  n/4
  58                 s += n.to_s + "/4"
  59           elsif ((a * 6) - (n = (a * 6 + 0.5).floor)).abs < 1e-4
  60             #  n/6
  61                 s += n.to_s + "/6"
  62           else
  63             s += sprintf("%.4f", a)
  64           end
  65         end
  66         ["x", "y", "z"].each_with_index { |t, i|
  67           a = sym[i, j]
  68           if a.abs < 1e-4
  69             next
  70           elsif (a - 1.0).abs < 1e-4
  71             s += (s == "" ? t : "+" + t)
  72           elsif (a + 1.0).abs < 1e-4
  73             s += "-" + t
  74           else
  75             if a > 0.0 && s != ""
  76                   s += "+"
  77                 elsif a < 0.0
  78                   s += "-"
  79                 end
  80                 s += sprintf("%.4f", a.abs) + t
  81           end
  82         }
  83         str += s
  84         if j < 2
  85           str += ", "
  86         end
  87   }
  88   str
  89 end
  90
  91 def string_to_symmetry(str)
  92   begin
  93     sary = str.split(/, */)
  94     raise if sary.count != 3
  95         a = []
  96         sary.each_with_index { |s, j|
  97           if s[0] != "-"
  98             s = "+" + s
  99           end
 100           while s.length > 0
 101             raise if s !~ /^([-+][.0-9\/]*)([xyzXYZ]?)/
 102                 sa = Regexp.last_match[1]
 103                 st = Regexp.last_match[2]
 104                 s = Regexp.last_match.post_match
 105             case st
 106                 when "x", "X"
 107                   i = 0
 108                 when "y", "Y"
 109                   i = 1
 110                 when "z", "Z"
 111                   i = 2
 112                 else
 113                   i = 3
 114                 end
 115                 raise if a[i * 3 + j] != nil
 116                 if sa == "-"
 117                   aa = -1
 118                 elsif sa == "+"
 119                   aa = 1
 120                 elsif sa.index("/")
 121                   sa0, sa1 = sa.split("/")
 122                   aa = sa0.to_f / sa1.to_f
 123                 else
 124                   aa = sa.to_f
 125                 end
 126                 a[i * 3 + j] = aa
 127           end
 128         }
 129         12.times { |i|
 130           a[i] ||= 0.0
 131         }
 132         return Transform.new(a)
 133   rescue
 134     raise "Cannot convert to symmetry operation: #{str}"
 135   end
 136 end
 137
 138 end
 139
 140 class Molecule
 141
 142 #  Export ortep to File fp
 143 #  If attr is a hash, it represents options for drawing.
 144 #  "atoms"=>[[group, type, color, rad], [group, type, color, rad], ...]
 145 #    group is an IntGroup, representing a group of atoms.
 146 #    type is an atom type: 0 = boundary, 1 = boundary + principal, 2 = 1 + axes, 3 = 2 + shades, 4 = fixed size
 147 #    color is 0..7 (0,1=black, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 7=yellow)
 148 #    rad is the radius of the atom (in Angstrom) whose type is fixed size
 149 #    If an atom appears in multiple entries, the later entry is valid.
 150 #  "bonds"=>[[group1, group2, type, color], [group1, group2, type, color], ...]
 151 #    group1 and group2 are IntGroups, reprensenting the atoms constituting the bonds.
 152 #    type is a bond type: 0 to 4 for bonds with no shades to 4 shades.
 153 #    color is 0..7 as in atoms.
 154 #    If a bond appears in multiple entries, the later entry is valid.
 155 def export_ortep(fp, attr = nil)
 156
 157   #  Create atom list
 158   hidden = atom_group { |ap| !is_atom_visible(ap.index) }
 159   hydrogen = self.show_hydrogens
 160   expanded = self.show_expanded
 161   atomlist = atom_group { |ap|
 162     (ap.element != "H" || hydrogen) &&
 163     (ap.symop == nil || expanded) &&
 164     (!hidden.include?(ap.index))
 165   }
 166
 167   #  Title
 168   fp.printf "%-78.78s\n", self.name + ": generated by Molby at " + Time.now.to_s
 169
 170   #  Cell parameters
 171   cp = self.cell
 172   if cp == nil
 173     cp = [1, 1, 1, 90, 90, 90]
 174   end
 175   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]
 176
 177   #  Symmetry operations
 178   syms = self.symmetries
 179   if syms == nil || syms.length == 0
 180    fp.print "1             0  1  0  0              0  0  1  0              0  0  0  1\n"
 181   else
 182     syms.each_with_index { |s, i|
 183       a = s.to_a
 184       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]
 185     }
 186   end
 187
 188   #  Atoms (all symmetry unique atoms regardless they are visible or not)
 189   n = 0
 190   aattr = (attr ? attr["atoms"] : nil)
 191   each_atom { |ap|
 192     break if ap.symop != nil
 193     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
 194         rad = -1.0
 195         if aattr
 196           aattr.reverse_each { |at|
 197             if at[0].include?(ap.index)
 198                   if at[1] == 4
 199                     rad = at[3].to_f
 200                     if rad == 0.0
 201                       rad = 0.1
 202                     end
 203                   end
 204                   break
 205                 end
 206           }
 207         end
 208     an = ap.aniso
 209     if an != nil
 210       eigval = ap.aniso_eigenvalues
 211       if eigval && (eigval[0] < 0.0 || eigval[1] < 0.0 || eigval[2] < 0.0)
 212           #  Non positive-definite anisotropic factor: fallback to isotropic
 213           an = nil
 214       end
 215     end
 216     if an != nil && rad < 0.0
 217       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
 218     else
 219           type = 7
 220           if rad < 0.0
 221         rad = ap.temp_factor
 222         rad = 1.2 if rad <= 0
 223             type = 6
 224           end
 225       fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", rad, 0.0, 0.0, 0.0, 0.0, 0.0, type
 226     end
 227     n += 1
 228   }
 229   natoms_tep = n
 230
 231   #  Special points to specify cartesian axes
 232   axis, angle = self.get_view_rotation
 233   tr = Transform.rotation(axis, -angle)
 234   org = self.get_view_center
 235   x = org + tr.column(0)
 236   y = org + tr.column(1)
 237   tr = self.cell_transform
 238   if tr
 239     tr = tr.inverse
 240   else
 241     tr = Transform.identity
 242   end
 243   org = tr * org
 244   x = tr * x
 245   y = tr * y
 246   fp.printf " CNTR                      %9.4f%9.4f%9.4f        0\n", org.x, org.y, org.z
 247   fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
 248   fp.printf " X                         %9.4f%9.4f%9.4f        0\n", x.x, x.y, x.z
 249   fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
 250   fp.printf " Y                         %9.4f%9.4f%9.4f        0\n", y.x, y.y, y.z
 251   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
 252
 253   #  Initialize
 254   fp.print  "      201\n"
 255
 256   #  Pen size
 257   fp.print  "      205        8\n"
 258
 259   #  Paper size, margin, viewing distance
 260   fp.print  "      301      6.6      6.6        0      0.0\n"
 261
 262   #  Sort atoms by symop and index
 263   acodes = Hash.new
 264   ig = IntGroup.new   #  Used for collecting contiguous atoms
 265   each_atom(atomlist) { |ap|
 266     idx, symcode = ap.to_adc
 267     acode = symcode * self.natoms + idx - 1
 268     acodes[acode] = ap.index
 269     ig.add(acode)
 270   }
 271   index2an = []       #  Index in Molecule to Index in ORTEP
 272   ig.each_with_index { |acode, i|
 273     index2an[acodes[acode]] = i + 1
 274   }
 275   i = 0
 276   adcs = []
 277   while (r = ig.range_at(i)) != nil
 278     s = r.first
 279     s = (s / self.natoms) + (s % self.natoms + 1) * 100000  #  Rebuild true ADC (atom index is in the upper digit)
 280     e = r.last
 281     e = (e / self.natoms) + (e % self.natoms + 1) * 100000
 282     if s < e
 283       adcs.push(s)
 284       adcs.push(-e)
 285     else
 286       adcs.push(s)
 287     end
 288     i += 1
 289   end
 290   k = 0
 291
 292   #  Atom list
 293   adcs.each_with_index { |a, i|
 294     if k == 0
 295       fp.print "      401"
 296     end
 297     fp.printf "%9d", a
 298     k += 1
 299     if i == adcs.length - 1 || k == 6 || (k == 5 && i < adcs.length - 2 && adcs[i + 2] < 0)
 300       fp.print "\n"
 301       k = 0
 302     end
 303   }
 304
 305   #  Axes
 306   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
 307 #  fp.print  "      502        1      0.0        2      0.0        3      0.0\n"
 308
 309   #  Autoscale
 310   fp.print  "      604                               1.538\n"
 311
 312   #  Explicit bonds
 313   bond_inst = Array.new(35) { [] }   #  Bonds for types 0 to 4 and colors 1 to 7
 314   battr = (attr ? attr["bonds"] : nil)
 315   bonds.each { |b|
 316     next if !atomlist.include?(b[0]) || !atomlist.include?(b[1])
 317         btype = bcol = nil
 318         if battr
 319           battr.reverse_each { |at|
 320             if (at[0].include?(b[0]) && at[1].include?(b[1])) || (at[0].include?(b[1]) && at[1].include?(b[0]))
 321                   btype = at[2] % 5
 322                   bcol = (at[3] != 0 ? at[3] - 1 : at[3]) % 7
 323                   break
 324                 end
 325           }
 326         end
 327         if btype == nil
 328           bcol = 0
 329       an1 = atoms[b[0]].atomic_number
 330       an2 = atoms[b[1]].atomic_number
 331       if an1 == 1 || an2 == 1
 332         btype = 0
 333       elsif an1 <= 8 && an2 <= 8
 334         btype = 2
 335       else
 336         btype = 4
 337       end
 338         end
 339     bond_inst[btype * 7 + bcol].push(b[0], b[1])
 340   }
 341
 342   #  Output bond specifications
 343   #  Avoid including too many ADCs in a single 811/821 instruction
 344   #  (Upper limit is 140. Here we divide at every 36 ADCs)
 345   output_bonds = lambda { |icode|
 346     bond_inst.each_with_index { |inst, ii|
 347       next if inst.length == 0
 348           btype = ii / 7 + 1
 349           if icode / 10 == 81   #  811 instructions
 350             fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
 351           end
 352       inst.each_with_index { |b, i|
 353         if i % 6 == 0
 354           fp.printf "  %d   %3s", (i >= inst.length - 6 || i % 36 == 30 ? 2 : 1), (i % 36 == 0 ? icode.to_s : "")
 355         end
 356         idx, scode = atoms[b].to_adc
 357         fp.printf "%9d", idx * 100000 + scode
 358         if i % 6 == 5 || i == inst.length - 1
 359           fp.print "\n"
 360           if i == inst.length - 1 || i % 36 == 35
 361             fp.printf "%21s%3d%12s%6.3f\n", "", btype, "", 0.05
 362           end
 363         end
 364       }
 365     }
 366   }
 367
 368   fp.print "  0  1001     0.02\n"  #  Activate hidden line removal
 369   output_bonds.call(821)
 370
 371   #  Atom types
 372   #  Atom type 0=714 (boundary), 1=712 (+principal), 2=716 (+axes), 3=711 (+shades)
 373   atom_inst = Array.new(28) { IntGroup.new }
 374   aattr = (attr ? attr["atoms"] : nil)
 375   atomlist.each { |i|
 376     atype = acol = nil
 377         if aattr
 378           aattr.reverse_each { |at|
 379             if at[0].include?(i)
 380                   atype = at[1] % 4
 381                   acol = (at[2] != 0 ? at[2] - 1 : at[2]) % 7
 382                   break
 383                 end
 384           }
 385         end
 386         if atype == nil
 387           acol = 0
 388       an1 = atoms[i].atomic_number
 389       if an1 == 1
 390         atype = 0
 391       elsif an1 <= 6
 392         atype = 1
 393       else
 394         atype = 3
 395       end
 396         end
 397     idx, scode = atoms[i].to_adc
 398     atom_inst[atype * 7 + acol].add(idx)
 399   }
 400   (atom_inst.count - 1).downto(0) { |ii|
 401     inst = atom_inst[ii]
 402         next if inst.count == 0
 403         fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
 404     i = 0
 405         atype = [714, 712, 716, 711][ii / 7]
 406     while (r = inst.range_at(i)) != nil
 407       fp.printf "  1   %3d\n", atype
 408       fp.printf "%27s%9d%9d\n", "", r.first, r.last
 409       i += 1
 410     end
 411   }
 412
 413   output_bonds.call(811)
 414
 415   #  Close plot
 416   fp.print "      202\n"
 417   fp.print "       -1\n"
 418
 419 end
 420
 421 def savetep(filename)
 422   if natoms == 0
 423         raise MolbyError, "cannot save ORTEP input; the molecule is empty"
 424   end
 425   fp = open(filename, "wb")
 426   export_ortep(fp)
 427   fp.close
 428   return true
 429 end
 430
 431 end
 432
 433 #  Best-fit planes
 434 #  Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
 435 #       T. Ito, Acta Cryst 1981, A37, 621-624
 436
 437 #  An object to describe the best-fit plane
 438 #
 439 #  Contains the plane coefficients and the constant (a, b, c, d for ax + by + cz + d = 0),
 440 #  the error matrix, and the metric tensor.
 441 #
 442 class Molby::Plane
 443   attr_accessor :molecule, :group, :coeff, :const, :error_matrix, :metric_tensor
 444   def initialize(mol, group, coeff, const, err, met)
 445     @molecule = mol
 446     @group = group
 447     @coeff = coeff
 448     @const = const
 449     @error_matrix = err
 450     @metric_tensor = met
 451     self
 452   end
 453   def coeffs
 454     [@coeff.x, @coeff.x, @coeff.z, @const]
 455   end
 456   def sigma
 457     [sqrt_safe(@error_matrix[0, 0]), sqrt_safe(@error_matrix[1, 1]), sqrt_safe(@error_matrix[2, 2]), sqrt_safe(@error_matrix[3, 3])]
 458   end
 459   def inspect
 460     s = sprintf("Molby::Plane[\n coeff, const = [[%f, %f, %f], %f],\n", @coeff.x, @coeff.y, @coeff.z, @const)
 461     s += sprintf(" sigma = [[%10.4e, %10.4e, %10.4e], %10.4e],\n", *self.sigma)
 462     (0..3).each { |i|
 463       s += (i == 0 ? " error_matrix = [" : "     ")
 464       (0..i).each { |j|
 465         s += sprintf("%12.6e%s", @error_matrix[j, i], (j == i ? (i == 3 ? "],\n" : ",\n") : ","))
 466       }
 467     }
 468     s += sprintf(" molecule = %s\n", @molecule.inspect)
 469     s += sprintf(" group = %s\n", @group.inspect)
 470     (0..3).each { |i|
 471       s += (i == 0 ? " metric_tensor = [" : "     ")
 472       (0..3).each { |j|
 473         s += sprintf("%12.6e%s", @metric_tensor[j, i], (j == 3 ? (i == 3 ? "]]\n" : ",\n") : ","))
 474       }
 475     }
 476     s
 477   end
 478   def distance(ap)
 479     if ap.is_a?(AtomRef)
 480       fr = ap.fract_r
 481       sig = ap.sigma
 482     else
 483       fr = Vector3D[*ap]
 484       sig = Vector3D[0, 0, 0]
 485     end
 486     d = fr.dot(@coeff) + @const
 487     sig1 = (@coeff.x * sig.x) ** 2 + (@coeff.y * sig.y) ** 2 + (@coeff.z * sig.z) ** 2
 488     sig2 = LAMatrix.multiply("t", fr, @error_matrix, fr)[0, 0]
 489     if ap.is_a?(AtomRef) && ap.molecule == @molecule && @group.include?(ap.index)
 490       #  The atom defines the plane
 491       sig0 = sig1 - sig2
 492       sig0 = 0.0 if sig0 < 0.0
 493     else
 494       sig0 = sig1 + sig2
 495     end
 496     return d, sqrt_safe(sig0)
 497   end
 498   def dihedral(plane)
 499     e1 = @error_matrix.submatrix(0, 0, 3, 3)
 500     e2 = plane.error_matrix.submatrix(0, 0, 3, 3)
 501     m = @metric_tensor.submatrix(0, 0, 3, 3)
 502     c = plane.coeff
 503     cos_t = plane.coeff.dot(m * @coeff)
 504     if cos_t > 1.0
 505       cos_t = 1.0
 506     elsif cos_t < -1.0
 507       cos_t = -1.0
 508     end
 509     t = acos(cos_t)
 510     sig_t = (m * e1).trace + (m * e2).trace
 511     if sig_t < t * t
 512       w = 1.0 / sin(t)
 513       sig_t = w * w * (c.dot(LAMatrix.multiply(m, e1, m, c)) + @coeff.dot(LAMatrix.multiply(m, e2, m, @coeff)))
 514     end
 515     t *= 180.0 / PI
 516     sig_t = sqrt_safe(sig_t) * 180.0 / PI
 517     return t, sig_t
 518   end
 519 end
 520
 521 class Molecule
 522
 523 #  Calculate best-fit plane for the given atoms
 524 #  Return value: a Molby::Plane object
 525
 526 def plane(group)
 527
 528   #  Number of atoms
 529   dim = group.length
 530
 531   #  Positional parameters and standard deviations
 532   x = []
 533   sig = []
 534   sig_min = 1e10
 535   each_atom(group) { |ap|
 536     x.push(ap.fract_r)
 537     sig.push(ap.sigma)
 538     if (s = ap.sigma_x) > 0.0 && s < sig_min
 539       sig_min = s
 540     end
 541     if (s = ap.sigma_y) > 0.0 && s < sig_min
 542       sig_min = s
 543     end
 544     if (s = ap.sigma_z) > 0.0 && s < sig_min
 545       sig_min = s
 546     end
 547   }
 548   if sig_min == 1e10
 549     sig_min = 1e-12
 550   end
 551   sig.each { |s|
 552     s.x = sig_min if s.x < sig_min
 553     s.y = sig_min if s.y < sig_min
 554     s.z = sig_min if s.z < sig_min
 555   }
 556
 557   #  The metric tensor of the reciprocal lattice
 558   #  g[j, i] = (ai*).dot(aj*), where ai* and aj* are the reciprocal axis vectors
 559   t = self.cell_transform
 560   if t.nil?
 561     t = Transform.identity
 562   end
 563   g2inv = LAMatrix[t]
 564   g2 = g2inv.inverse
 565   g2[3, 3] = 0.0
 566   g2inv[3, 3] = 0.0
 567   g = LAMatrix.multiply("t", g2, g2)
 568
 569   #  The variance-covariance matrices of the atomic parameters
 570   #  mm[k][n] is a 3x3 matrix describing the correlation between the atoms k and n,
 571   #  and its components are defined as: sigma_k[i] * sigma_n[j] * corr[k, i, n, j],
 572   #  where corr(k, i, n, j) is the correlation coefficients between the atomic parameters
 573   #  k[i] and n[j].
 574   mm = Array.new(dim) { Array.new(dim) }
 575   zero = LAMatrix.zero(3, 3)
 576   dim.times { |k|
 577     dim.times { |n|
 578       mkn = LAMatrix.new(3, 3)
 579       if k == n
 580         3.times { |i|
 581           3.times { |j|
 582             if i == j
 583               mkn[j, i] = sig[k][i] * sig[n][j]
 584             else
 585               #  Inter-coordinate correlation should be implemented here
 586             end
 587           }
 588         }
 589       else
 590         #  Inter-atomic correlation should be implemented here
 591       end
 592       mm[k][n] = (mkn == zero ? zero : mkn)
 593     }
 594   }
 595
 596   #  The variance-covariance matrix of the atom-plance distances
 597   #  m[j, i] = v.transpose * mm[i][j] * v, where v is the plane coefficient vector
 598   #  The inverse of m is the weight matrix
 599   m = LAMatrix.new(dim, dim)
 600
 601   #  The matrix representation of the atomic coordinates
 602   #  y[j, i] = x[i][j] (for j = 0..2), -1 (for j = 3)
 603   #  y * LAMatrix[a, b, c, d] gives the atom-plane distances for each atom
 604   y = LAMatrix.new(4, dim)
 605   dim.times { |i|
 606     y[0, i] = x[i].x
 607     y[1, i] = x[i].y
 608     y[2, i] = x[i].z
 609     y[3, i] = 1.0
 610   }
 611
 612   #  The coefficients to be determined
 613   n0 = LAMatrix[1, 1, 1, 0]
 614   v = LAMatrix[1, 1, 1]     #  The coefficient part
 615
 616   iter = 0
 617   while iter < 20
 618
 619     iter += 1
 620
 621     #  Set zero to the "constant" part, and normalize the "coefficient" part
 622     n0[0, 3] = 0.0
 623     n0 = g2 * n0
 624     n0.multiply!(1.0 / n0.fnorm)
 625     n0 = g2inv * n0
 626     3.times { |i| v[0, i] = n0[0, i] }
 627
 628     #  Build the variance-covariance matrix
 629     dim.times { |i|
 630       dim.times { |j|
 631         m[j, i] = LAMatrix.multiply("t", v, mm[i][j], v)[0, 0]
 632       }
 633     }
 634     c = LAMatrix.multiply("t", y, "i", m, y)
 635
 636     #  Invert c: only the inverse is used in the following, so c is inversed destructively
 637     cinv = c.inverse!
 638
 639     if iter == 1
 640
 641       #  Determine the tentative solution, which is given by the eigenvector of cinv * g
 642       #  for the largest eigenvalue
 643       evals, evecs = (cinv * g).eigenvalues
 644       4.times { |i| n0[0, i] = evecs[3, i] }
 645
 646     else
 647
 648       #  Convert the coefficient vector to the reciprocal space
 649       h = g * n0
 650
 651       #  Determine multiplier
 652       #  In this implementation, the sign of delta-n is opposite from that used in
 653       #  the reference
 654       lam = 1.0 / (LAMatrix.multiply("t", h, cinv, h)[0, 0])
 655
 656       #  Solve the linearized equation
 657       #  (Is the equation 21 in the reference really correct? Shouldn't it read
 658       #   B = 1 - lambda * C.inverse * H* ? )
 659       b = LAMatrix.multiply(lam, cinv, g)
 660       b.sub!(LAMatrix.identity(4))
 661
 662       dn = b * n0
 663       n0 += dn
 664
 665       break if dn[0, 0] ** 2 + dn[0, 1] ** 2 + dn[0, 2] ** 2 < 1e-9
 666
 667     end
 668   end
 669
 670   #  Error matrix = b * cinv * b.transpose
 671   em = LAMatrix.multiply(b, cinv, "t", b)
 672   coeff = Vector3D[n0[0, 0], n0[0, 1], n0[0, 2]]
 673   const = n0[0, 3]
 674
 675   return Molby::Plane.new(self, group, coeff, const, em, g)
 676
 677 end
 678
 679 def cmd_plane
 680   plane_settings = @plane_settings || Hash.new
 681   mol = self
 682   mol.open_auxiliary_window("Best-Fit Planes", :has_close_box=>true) {
 683     refresh_proc = lambda { |it|
 684       n = it[:tag][/\d/].to_i
 685       g = plane_settings["group#{n}"]
 686       if g
 687         str = g.inspect.sub!("IntGroup[", "").sub!("]", "")
 688         set_value("group#{n}", str)
 689         if n == 1 || n == 2
 690           p = mol.plane(g) rescue p = nil
 691           plane_settings["plane#{n}"] = p
 692           if p
 693             coeff = p.coeff
 694             const = p.const
 695             sig = p.sigma
 696             aps = (n == 1 ? "" : "'")
 697             str = sprintf("a%s = %f(%f)\nb%s = %f(%f)\nc%s = %f(%f)\nd%s = %f(%f)",
 698                           aps, coeff.x, sig[0],
 699                           aps, coeff.y, sig[1],
 700                           aps, coeff.z, sig[2],
 701                           aps, const, sig[3])
 702             set_value("result#{n}", str)
 703           else
 704             set_value("result#{n}", "")
 705           end
 706           p1 = plane_settings["plane1"]
 707           p2 = plane_settings["plane2"]
 708           if p1 && p2
 709             t, sig = p1.dihedral(p2)
 710             str = sprintf("%f(%f)", t, sig)
 711             set_value("dihedral", str)
 712           else
 713             set_value("dihedral", "")
 714           end
 715         else
 716           p = plane_settings["plane1"]
 717           if p
 718             str = ""
 719             mol.each_atom(g) { |ap|
 720               d, sig = p.distance(ap)
 721               str += sprintf("%s  %f(%f)\n", ap.name, d, sig)
 722             }
 723             str.chomp!
 724           else
 725             str = ""
 726           end
 727           set_value("result#{n}", str)
 728         end
 729       else
 730         set_value("group#{n}", "")
 731         set_value("result#{n}", "")
 732       end
 733     }
 734     set_proc = lambda { |it|
 735       n = it[:tag][/\d/].to_i
 736       sel = mol.selection
 737       if sel.count > 0
 738         str = sel.inspect.sub!("IntGroup[", "").sub!("]", "")
 739         set_value("group#{n}", str)
 740         plane_settings["group#{n}"] = sel
 741       else
 742         plane_settings["group#{n}"] = nil
 743       end
 744       refresh_proc.call(it)
 745     }
 746     text_proc = lambda { |it|
 747       n = it[:tag][/\d/].to_i
 748       str = it[:value].gsub(/[^-.,0-9]/, "")  #  Remove unsane characters
 749       g = eval("IntGroup[#{str}]") rescue g = nil
 750       plane_settings["group#{n}"] = g
 751       refresh_proc.call(it)
 752     }
 753     layout(3,
 754       item(:text, :title=>"Plane 1 (ax + by + cz + d = 0)"),
 755       -1, -1,
 756       item(:text, :title=>"Atoms"),
 757       item(:textfield, :width=>240, :height=>32, :tag=>"group1", :action=>text_proc),
 758       item(:button, :title=>"Set Current Selection", :tag=>"button1", :action=>set_proc),
 759       item(:text, :title=>"Results"),
 760       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result1"),
 761       item(:button, :title=>"Recalculate", :tag=>"refresh1", :action=>refresh_proc),
 762       item(:line),
 763       -1, -1,
 764       item(:text, :title=>"Plane 2 (a'x + b'y + c'z + d' = 0)"),
 765       -1, -1,
 766       item(:text, :title=>"Atoms"),
 767       item(:textfield, :width=>240, :height=>32, :tag=>"group2", :action=>text_proc),
 768       item(:button, :title=>"Set Current Selection", :tag=>"button2", :action=>set_proc),
 769       item(:text, :title=>"Results"),
 770       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result2"),
 771       item(:button, :title=>"Recalculate", :tag=>"refresh2", :action=>refresh_proc),
 772       item(:text, :title=>"Dihedral angle with Plane 1"), -1, -1,
 773       -1,
 774       item(:textfield, :width=>240, :height=>20, :tag=>"dihedral"), -1,
 775       item(:line),
 776       -1, -1,
 777       item(:text, :title=>"Distance from Plane 1"), -1, -1,
 778       item(:text, :title=>"Atoms"),
 779       item(:textfield, :width=>240, :height=>32, :tag=>"group3", :action=>text_proc),
 780       item(:button, :title=>"Set Current Selection", :tag=>"button3", :action=>set_proc),
 781       item(:text, :title=>"Results"),
 782       item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result3"),
 783       item(:button, :title=>"Recalculate", :tag=>"refresh3", :action=>refresh_proc)
 784     )
 785     refresh_proc.call(item_with_tag("refresh1"))
 786     refresh_proc.call(item_with_tag("refresh2"))
 787     refresh_proc.call(item_with_tag("refresh3"))
 788         show
 789   }
 790   @plane_settings = plane_settings
 791 end
 792
 793 #  Calculate bond length and angles with standard deviations
 794 #  args can be a single IntGroup (calculate bonds and angles including those atoms)
 795 #  or arrays of 2 or 3 integers (explicitly specifying bonds and angles)
 796 def bond_angle_with_sigma(*args)
 797
 798   if args.length >= 2 || (args.length == 1 && !args[0].is_a?(IntGroup))
 799     #  Bonds and angles are explicitly specified
 800     bonds = []
 801     angles = []
 802     args.each { |arg|
 803       if arg.length == 2 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer)
 804         bonds.push(arg)
 805       elsif arg.length == 3 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer) && arg[2].is_a?(Integer)
 806         angles.push(arg)
 807       else
 808         raise MolbyError, "Invalid argument #{arg.inspect}"
 809       end
 810     }
 811   else
 812     if args.length == 0
 813           g = nil
 814         else
 815           g = args[0]
 816         end
 817         bonds = self.bonds.select { |b|
 818           (g == nil || (g.include?(b[0]) && g.include?(b[1]))) &&
 819             (atoms[b[0]].symop == nil || atoms[b[1]].symop == nil)
 820         }
 821         angles = self.angles.select { |ang|
 822           (g == nil || (g.include?(ang[0]) && g.include?(ang[1]) && g.include?(ang[2]))) &&
 823             (atoms[ang[0]].symop == nil || atoms[ang[1]].symop == nil || atoms[ang[2]].symop == nil)
 824         }
 825   end
 826
 827   #  A list of interatomic distance, its differential coefficients, and other
 828   #  useful quantities.
 829   #  The hash tag is a list [i, j] (i and j are the atom indices) or [i, j, k]
 830   #  (i, j, k are the atom indices, and r(ijk) is the distance between the atom i
 831   #  and the center point between the atoms j and k).
 832   #  The value is a list of following quantities:
 833   #    index 0: rij or r(ijk)
 834   #    index 1-9: d(rij)/d(xij), d(rij)/d(yij), d(rij)/d(zij),
 835   #      d(rij)/da, d(rij)/db, d(rij)/dc, d(rij)/d(alpha), d(rij)/d(beta), d(rij)/d(gamma)
 836   #    index 10: the list of the "base atom"
 837   #    index 11: the list of the transform matrices
 838   dlist = Hash.new
 839
 840   #  A list of fractional coordinates and sigmas
 841   fract = []
 842   sigma = []
 843   each_atom { |ap|
 844     fract.push(ap.fract_r)
 845     sigma.push(ap.sigma)
 846   }
 847
 848   #  A list of base atoms (for symmetry-related atoms) and transform matrices
 849   bases = []
 850   trans = []
 851   symcode = []
 852   trans_i = Transform.identity
 853   each_atom { |ap|
 854     sym = ap.symop
 855     bases.push(sym ? sym[4] : ap.index)
 856         if sym
 857           tr = transform_for_symop(sym).transpose
 858       tr[3, 0] = tr[3, 1] = tr[3, 2] = 0.0
 859         else
 860           tr = trans_i
 861         end
 862     trans.push(tr)
 863         symcode.push(sym ? sprintf("%d_%d%d%d", sym[0] + 1, sym[1] + 5, sym[2] + 5, sym[3] + 5) : ".")
 864   }
 865
 866   #  Unit cell parameter
 867   cell = self.cell
 868   if cell == nil
 869     cell = [1, 1, 1, 90, 90, 90, 0, 0, 0, 0, 0, 0]
 870   elsif cell.length == 6
 871     cell.push(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)  #  No sigma information
 872   end
 873   # $a, $b, $c, $alpha, $beta, $gamma, $sig_a, $sig_b, $sig_c, $sig_alpha, $sig_beta, $sig_gamma = self.cell
 874   cos_a = cos(cell[3] * PI / 180.0)
 875   cos_b = cos(cell[4] * PI / 180.0)
 876   cos_c = cos(cell[5] * PI / 180.0)
 877   abc = cell[0] * cell[1] * cos_c
 878   bca = cell[1] * cell[2] * cos_a
 879   cab = cell[2] * cell[0] * cos_b
 880   aa = cell[0] * cell[0]
 881   bb = cell[1] * cell[1]
 882   cc = cell[2] * cell[2]
 883
 884   get_dlist = lambda { |_i, _j, _k|
 885     if _k != nil
 886       if _j > _k
 887         _j, _k = _k, _j
 888       end
 889       _p = dlist[[_i, _j, _k]]
 890       return _p if _p != nil
 891       _p = (dlist[[_i, _j, _k]] = [])
 892       _vij = fract[_i] - (fract[_j] + fract[_k]) * 0.5
 893       _p[10] = [bases[_i], bases[_j], bases[_k]]
 894       _p[11] = [trans[_i], trans[_j], trans[_k]]
 895     else
 896       if _i > _j
 897         _i, _j = _j, _i
 898       end
 899       _p = dlist[[_i, _j]]
 900           return _p if _p != nil
 901       _p = (dlist[[_i, _j]] = [])
 902       _vij = fract[_i] - fract[_j]
 903       _p[10] = [bases[_i], bases[_j]]
 904       _p[11] = [trans[_i], trans[_j]]
 905     end
 906     _xij = _vij.x
 907     _yij = _vij.y
 908     _zij = _vij.z
 909     _dij = sqrt_safe(aa * _xij * _xij + bb * _yij * _yij + cc * _zij * _zij + 2 * bca * _yij * _zij + 2 * cab * _zij * _xij + 2 * abc * _xij * _yij)
 910     _p[0] = _dij
 911     _p[1] = (aa * _xij + abc * _yij + cab * _zij) / _dij
 912     _p[2] = (bb * _yij + bca * _zij + abc * _xij) / _dij
 913     _p[3] = (cc * _zij + cab * _xij + bca * _yij) / _dij
 914     _p[4] = (cell[0] * _xij * _xij + cell[2] * cos_b * _zij * _xij + cell[1] * cos_c * _xij * _yij) / _dij
 915     _p[5] = (cell[1] * _yij * _yij + cell[0] * cos_c * _xij * _yij + cell[2] * cos_a * _yij * _zij) / _dij
 916     _p[6] = (cell[2] * _zij * _zij + cell[1] * cos_a * _yij * _zij + cell[0] * cos_b * _zij * _xij) / _dij
 917     _p[7] = (-cell[1] * cell[2] * sin(cell[3] * PI / 180.0) * _yij * _zij) * (PI / 180.0) / _dij
 918     _p[8] = (-cell[2] * cell[0] * sin(cell[4] * PI / 180.0) * _zij * _xij) * (PI / 180.0) / _dij
 919     _p[9] = (-cell[0] * cell[1] * sin(cell[5] * PI / 180.0) * _xij * _yij) * (PI / 180.0) / _dij
 920     return _p
 921   }
 922
 923   diff_by_rn = lambda { |_dl, _n, _ijk|
 924     #  dl = dlist(i, j)
 925     #  return value: Vector3D[ d(rij)/d(xn), d(rij)/d(yn), d(rij)/d(zn) ]
 926     #  If ijk is true, then dl is dlist(i, j, k)
 927     _dv = Vector3D[_dl[1], _dl[2], _dl[3]]
 928     _c = Vector3D[0, 0, 0]
 929     if _dl[10][0] == _n
 930       _c += _dl[11][0] * _dv
 931     end
 932     if _ijk
 933       if _dl[10][1] == _n
 934         _c -= _dl[11][1] * _dv * 0.5
 935       end
 936       if _dl[10][2] == _n
 937         _c -= _dl[11][2] * _dv * 0.5
 938       end
 939     else
 940       if _dl[10][1] == _n
 941         _c -= _dl[11][1] * _dv
 942       end
 943     end
 944         return _c
 945   }
 946
 947   notate_with_sigma = lambda { |_val, _sig|
 948     if _sig == 0.0
 949       return sprintf "%.4f", _val
 950     end
 951     _lg = log(_sig.abs / 1.95) / log(10.0)
 952     _n = -(_lg.floor)
 953     _val2 = (_val * (10 ** _n) + 0.5).floor * (0.1 ** _n)
 954     return sprintf "%.#{_n}f(%d)", _val2, (_sig * (10 ** _n) + 0.5).floor
 955   }
 956
 957   results = []
 958
 959   c = Vector3D[0, 0, 0]
 960   bonds.each { |b|
 961     i = b[0]
 962     j = b[1]
 963     if i > j
 964       i, j = j, i
 965     end
 966     p = get_dlist.call(i, j, nil)
 967     sig = 0.0
 968     p[10].uniq.each { |k|
 969       s = sigma[k]
 970       next unless s
 971       c = diff_by_rn.call(p, k, false)
 972       #  Test
 973       if nil
 974         apk = atoms[k]
 975         r = apk.fract_r
 976         d0 = calc_bond(i, j)
 977         apk.fract_r = r + Vector3D[0.00001, 0, 0]
 978         amend_by_symmetry(p[10])
 979         d1 = calc_bond(i, j)
 980         apk.fract_r = r + Vector3D[0, 0.00001, 0]
 981         amend_by_symmetry(p[10])
 982         d2 = calc_bond(i, j)
 983         apk.fract_r = r + Vector3D[0, 0, 0.00001]
 984         amend_by_symmetry(p[10])
 985         d3 = calc_bond(i, j)
 986         apk.fract_r = r
 987         amend_by_symmetry(p[10])
 988         printf " dr/dv[%d] cal %10.5g %10.5g %10.5g\n", k, c[0], c[1], c[2]
 989         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
 990       end
 991       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]
 992     }
 993     6.times { |n|
 994       sig += (p[4 + n] * cell[6 + n]) ** 2
 995     }
 996     sig = sqrt_safe(sig)
 997     results.push([[i, j], notate_with_sigma.call(p[0], sig), symcode[i], symcode[j], nil])
 998   }
 999
1000   angles.each { |ang|
1001     i = ang[0]
1002     j = ang[1]
1003     k = ang[2]
1004     p0 = get_dlist.call(i, j, nil)
1005     p1 = get_dlist.call(j, k, nil)
1006     p2 = get_dlist.call(i, k, nil)
1007     p3 = get_dlist.call(j, i, k)
1008     t123 = acos_safe((p0[0] ** 2 + p1[0] ** 2 - p2[0] ** 2) / (2 * p0[0] * p1[0]))
1009     t124 = acos_safe((p0[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p0[0] * p3[0]))
1010     t324 = acos_safe((p1[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p1[0] * p3[0]))
1011     dtdr12 = -(p0[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p0[0] ** 2) * sin(t124))
1012     dtdr23 = -(p1[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p1[0] ** 2) * sin(t324))
1013     dtdr13 = p2[0] / (sin(t124) * 4 * p0[0] * p3[0]) + p2[0] / (sin(t324) * 4 * p1[0] * p3[0])
1014     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))
1015     pp = (p0[10] + p1[10] + p2[10] + p3[10]).uniq
1016     sig = 0.0
1017     pp.each { |n|
1018       s = sigma[n]
1019       next unless s
1020       c = Vector3D[0, 0, 0]
1021       c += diff_by_rn.call(p0, n, false) * (dtdr12 * 180.0 / PI)
1022       c += diff_by_rn.call(p1, n, false) * (dtdr23 * 180.0 / PI)
1023       c += diff_by_rn.call(p2, n, false) * (dtdr13 * 180.0 / PI)
1024       c += diff_by_rn.call(p3, n, true) * (dtdr24 * 180.0 / PI)
1025       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]
1026     }
1027     dd = dtdr12 * p0[4] + dtdr23 * p1[4] + dtdr13 * p2[4] + dtdr24 * p3[4]
1028     sig += dd * dd * cell[6] * cell[6]
1029     dd = dtdr12 * p0[5] + dtdr23 * p1[5] + dtdr13 * p2[5] + dtdr24 * p3[5]
1030     sig += dd * dd * cell[7] * cell[7]
1031     dd = dtdr12 * p0[6] + dtdr23 * p1[6] + dtdr13 * p2[6] + dtdr24 * p3[6]
1032     sig += dd * dd * cell[8] * cell[8]
1033     dd = dtdr12 * p0[7] + dtdr23 * p1[7] + dtdr13 * p2[7] + dtdr24 * p3[7]
1034     sig += dd * dd * cell[9] * cell[9]
1035     dd = dtdr12 * p0[8] + dtdr23 * p1[8] + dtdr13 * p2[8] + dtdr24 * p3[8]
1036     sig += dd * dd * cell[10] * cell[10]
1037     dd = dtdr12 * p0[9] + dtdr23 * p1[9] + dtdr13 * p2[9] + dtdr24 * p3[9]
1038     sig += dd * dd * cell[11] * cell[11]
1039     sig = sqrt_safe(sig)
1040     results.push([[i, j, k], notate_with_sigma.call(t123*180/PI, sig), symcode[i], symcode[j], symcode[k]])
1041   }
1042   results
1043 end
1044
1045 def cmd_bond_angle_with_sigma
1046   mol = self
1047   mol.open_auxiliary_window("Bond & Angle with Sigma", :has_close_box=>true) {
1048     values = []
1049         clicked = []
1050         sel = mol.selection
1051         on_add_bond = lambda { |it|
1052           values.push([nil, nil, "-", nil, nil, nil, "-", 0])
1053           item_with_tag("table")[:refresh] = true
1054           item_with_tag("table")[:selection] = IntGroup[values.count - 1]
1055           item_with_tag("dump")[:enabled] = true
1056         }
1057         on_add_angle = lambda { |it|
1058           values.push([nil, nil, nil, nil, nil, nil, nil, 0])
1059           item_with_tag("table")[:refresh] = true
1060           item_with_tag("table")[:selection] = IntGroup[values.count - 1]
1061           item_with_tag("dump")[:enabled] = true
1062         }
1063         on_get_value = lambda { |it, row, col|
1064           v = values[row][col]
1065           if col < 3
1066             if v.is_a?(Integer)
1067               v = mol.atoms[v].name
1068                 elsif v == nil
1069                   v = "(select)"
1070                 end
1071           elsif (col >= 4 && col <= 6) && (vv = values[row][col - 4]).is_a?(Integer)
1072             if (s = mol.atoms[vv].symop) != nil
1073                   v = sprintf("%d_%d%d%d", s[0] + 1, s[1] + 5, s[2] + 5, s[3] + 5)
1074                 else
1075                   v = "."
1076                 end
1077           end
1078           return v
1079         }
1080         on_selection_changed = lambda { |it|
1081           s = it[:selection]
1082           if s && s.count > 0
1083             set_attr("delete", :enabled=>true)
1084             set_attr("dump", :enabled=>true)
1085           else
1086             set_attr("delete", :enabled=>false)
1087                 set_attr("dump", :enabled=>false)
1088           end
1089         }
1090         on_delete = lambda { |it|
1091           s = attr("table", :selection)
1092           if s
1093             s.reverse_each { |idx|
1094                   values.delete_at(idx)
1095                 }
1096           end
1097           tbl = item_with_tag("table")
1098           tbl[:selection] = IntGroup[]
1099           tbl[:refresh] = true
1100           on_selection_changed.call(tbl)
1101         }
1102         atom_name = lambda { |ap|
1103           (ap.molecule.nresidues >= 2 ? "#{ap.res_seq}:" : "") + ap.name
1104         }
1105     layout(1,
1106           item(:table, :width=>500, :height=>300, :tag=>"table",
1107             :columns=>[["atom1", 60], ["atom2", 60], ["atom3", 60],
1108                   ["value(sigma)", 80], ["symop1", 80], ["symop2", 80], ["symop3", 80]],
1109                 :on_count=> lambda { |it| values.count },
1110                 :on_get_value=>on_get_value,
1111                 :on_selection_changed=>on_selection_changed),
1112           item(:text, :title=>"(1) Hit 'New Bond' or 'New Angle', and (2) click on the atoms to measure in the main window.",
1113             :font=>[10]),
1114           layout(3,
1115             item(:button, :title=>"New Bond", :width=>80, :font=>[10], :action=>on_add_bond),
1116                 item(:button, :title=>"New Angle", :width=>80, :font=>[10], :action=>on_add_angle),
1117                 item(:button, :title=>"Delete", :width=>80, :font=>[10], :action=>on_delete,
1118                   :tag=>"delete", :enabled=>false),
1119                 :padding=>5, :margin=>0),
1120           layout(2,
1121             item(:view, :width=>480, :height=>1), -1,
1122             item(:button, :title=>"Export to Clipboard", :tag=>"dump",
1123                   :action=>lambda { |item|
1124                     s = ""
1125                         sel = attr("table", :selection)
1126                         return if sel == nil || sel.count == 0
1127                         sel.each { |j|
1128                           ss = ""
1129                           7.times { |i|
1130                             sss = on_get_value.call(item_with_tag("table"), j, i).to_s
1131                                 ss += (sss == "-" ? "" : sss) + (i == 6 ? "\n" : "  ")
1132                           }
1133                           s += ss
1134                         }
1135                         export_to_clipboard(s)
1136                   },
1137                   :enabled=>false),
1138                 [item(:button, :title=>"Close", :action=>lambda { |item| hide }), {:align=>:right}])
1139         )
1140         @on_document_modified = lambda { |*d|
1141           newsel = mol.selection - sel
1142           sel = mol.selection
1143           row = (item_with_tag("table")[:selection] || [nil])[-1]
1144           return unless row
1145           if newsel.count == 1
1146             #  Add new atom to the selected row
1147                 val = newsel[0]
1148                 lim = (values[row][2] == "-" ? 2 : 3)
1149                 i = 0
1150                 while i < lim
1151                   if values[row][i] == nil
1152                     values[row][i] = val
1153                         break
1154                   end
1155                   i += 1
1156                 end
1157                 if i == lim
1158                   (lim - 1).times { |j|
1159                     values[row][j] = values[row][j + 1]
1160                   }
1161                   values[row][lim - 1] = val
1162                 end
1163                 if i < lim - 1
1164                   val1 = nil
1165                 else
1166                   val1 = mol.bond_angle_with_sigma(values[row][0...lim])[0][1]
1167                 end
1168                 values[row][3] = val1
1169                 item_with_tag("table")[:refresh] = true
1170           end
1171         }
1172 #   listen(mol, "documentModified", on_document_modified)
1173 #       listen(mol, "documentWillClose", lambda { |*d| hide } )
1174         @on_document_modified.call
1175         show
1176   }
1177 end
1178
1179 def find_expanded_atom(base, symop)
1180   return base unless symop
1181   symopb = symop + [base]
1182   self.each_atom { |ap|
1183     if ap.symop == symopb
1184       return ap.index
1185     end
1186   }
1187   nil
1188 end
1189
1190 def complete_by_symmetry
1191   if self.box == nil
1192     raise "Unit cell should be given"
1193   end
1194   verbose = false
1195   avec, bvec, cvec = self.box
1196   syms = []
1197   self.nsymmetries.times { |n|
1198     syms.push(transform_for_symop([n, 0, 0, 0], true))
1199   }
1200   frags = []
1201   self.each_fragment { |f|
1202     frags.push(f)
1203   }
1204   close_pairs = []
1205   self.each_atom { |ap|
1206     #  Find if ap is contained in expansion
1207         next if ap.symop != nil  #  Already expanded
1208     rad = Parameter.builtin.elements[ap.atomic_number].radius
1209     syms.each_with_index { |sym, sym_idx|
1210       27.times { |n|
1211         dx = n % 3 - 1
1212         dy = (n / 3) % 3 - 1
1213         dz = (n / 9) % 3 - 1
1214         next if dx == 0 && dy == 0 && dz == 0 && sym_idx == 0
1215         symop = [sym_idx, dx, dy, dz]
1216                 next if self.find_expanded_atom(ap.index, symop)
1217         r = sym * ap.r + avec * dx + bvec * dy + cvec * dz
1218         close_atoms = self.find_close_atoms(r, 1.2, rad)
1219         if close_atoms.length > 0
1220           #  ap * [sym, dx, dy, dz] is included in the expansion
1221           close_atoms.each { |ca|
1222             next if ca > ap.index
1223             i1 = frags.index { |f| f.include?(ca) }
1224             i2 = frags.index { |f| f.include?(ap.index) }
1225             pp = close_pairs.find { |p1| p1[0] == i1 && p1[1] == i2 && p1[2] == symop }
1226                         if pp == nil
1227                           pp = [i1, i2, symop, [], [], [], []]
1228                           close_pairs.push(pp)
1229                     end
1230             if (r - atoms[ca].r).length2 > 0.0001
1231                           #  Normal case (bond between ca and ap)
1232                           pp[3].push(ca)
1233                           pp[4].push(ap.index)
1234                         else
1235               #  Special position
1236               pp[5].push(ca)
1237               pp[6].push(ap.index)
1238             end
1239           }
1240         end
1241       }
1242     }
1243   }
1244   puts "close_pairs = #{close_pairs}" if verbose
1245   expand_pairs = lambda { |i1, symop, depth|
1246     #  Find expanded fragment that is close to fragment [i1, symop]
1247     next [] if depth > 16
1248     retval = []
1249     close_pairs.each { |pp|
1250       next if i1 != nil && pp[0] != i1
1251           next if pp[3].count == 0   #  No real expansion
1252       #  Multiply two symops
1253       if symop
1254         symop2 = symop_for_transform(transform_for_symop(pp[2]) * transform_for_symop(symop))
1255         puts "multiply two symops: #{pp[2].inspect} * #{symop.inspect} -> #{symop2.inspect}" if verbose
1256       else
1257         symop2 = pp[2]
1258       end
1259       if symop2[0] == 0
1260         #  Translation only
1261         if symop2[1] == 0 && symop2[2] == 0 && symop2[3] == 0
1262                   #  No expansion, but pp[3] and pp[4] are still needed to make bonds
1263                   symop2 = nil
1264             end
1265                 #  Expand the atoms only at the marginal
1266             retval.push([nil, symop, symop2, pp[3], pp[4], pp[5], pp[6]])
1267       else
1268         #  Expand fragment
1269         retval.push([pp[1], symop, symop2, pp[3], pp[4], pp[5], pp[6]])
1270         retval.concat(expand_pairs.call(pp[1], symop2, depth + 1))
1271       end
1272     }
1273     next retval
1274   }
1275   ex_pairs = expand_pairs.call(nil, nil, 0)
1276   puts "ex_pairs = #{ex_pairs}" if verbose
1277
1278   #  Expand fragments
1279   duplicates = []
1280   ex_pairs.each { |ex|
1281     #  ex[0]: fragment to expand, ex[1], ex[2]: symop,
1282     #  ex[3], ex[4]: atoms to make bonds, ex[5], ex[6]: duplicate atoms
1283     n1 = self.natoms
1284     if ex[0] == nil
1285       f = IntGroup[ex[4]]
1286     else
1287       f = frags[ex[0]]
1288     end
1289         if ex[2] != nil
1290       self.expand_by_symmetry(f, ex[2][0], ex[2][1], ex[2][2], ex[2][3], true)
1291       puts "expand(#{f}, #{ex[2]}) -> atoms[#{n1}..#{self.natoms - n1}]" if verbose
1292         end
1293     ex[3].each_with_index { |x, i|
1294       x1 = self.find_expanded_atom(x, ex[1])
1295           if ex[2] != nil
1296         x2 = f.index(ex[4][i]) + n1  #  New index of the expanded atom
1297           else
1298             x2 = ex[4][i]   #  Base atom
1299           end
1300       create_bond(x1, x2)
1301       puts "create_bond(#{x1}, #{x2})" if verbose
1302     }
1303     #  Register duplicate atoms
1304         if ex[2] != nil
1305       ex[5].each_with_index { |x, i|
1306         x1 = self.find_expanded_atom(x, ex[1])
1307         x2 = f.index(ex[6][i]) + n1
1308         duplicates.each { |d|
1309           if d.include?(x1)
1310             d.push(x2)
1311             x2 = nil
1312             break
1313           end
1314         }
1315         if x2
1316           duplicates.push([x1, x2])
1317         end
1318       }
1319     end
1320   }
1321   puts "duplicates = #{duplicates}" if verbose
1322
1323   #  Remove duplicate atoms
1324   rem = []
1325   duplicates.each { |d|
1326     d.each_with_index { |dn, i|
1327       next if i == 0
1328       #  Remake bonds
1329       self.atoms[dn].connects.each { |nn|
1330         create_bond(d[0], nn)
1331         puts "create_bond(#{d[0]}, #{nn})" if verbose
1332       }
1333       rem.push(dn)
1334     }
1335   }
1336   remove(rem)
1337   puts "remove(#{rem})" if verbose
1338
1339 end
1340
1341 def create_packing_diagram
1342   if self.box == nil
1343     error_message_box "Unit cell is not defined."
1344         return
1345   end
1346   expansion_box = (@expansion_box ||= [0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
1347   h = Dialog.run("Create Packing Diagram") {
1348     layout(4,
1349           item(:text, :title=>"Specify the expansion limits for each axis:"),
1350           -1, -1, -1,
1351           item(:text, :title=>"x"),
1352           item(:textfield, :width=>80, :tag=>"xmin", :value=>sprintf("%.1f", expansion_box[0].to_f)),
1353           item(:text, :title=>"..."),
1354           item(:textfield, :width=>80, :tag=>"xmax", :value=>sprintf("%.1f", expansion_box[1].to_f)),
1355           item(:text, :title=>"y"),
1356           item(:textfield, :width=>80, :tag=>"ymin", :value=>sprintf("%.1f", expansion_box[2].to_f)),
1357           item(:text, :title=>"..."),
1358           item(:textfield, :width=>80, :tag=>"ymax", :value=>sprintf("%.1f", expansion_box[3].to_f)),
1359           item(:text, :title=>"z"),
1360           item(:textfield, :width=>80, :tag=>"zmin", :value=>sprintf("%.1f", expansion_box[4].to_f)),
1361           item(:text, :title=>"..."),
1362           item(:textfield, :width=>80, :tag=>"zmax", :value=>sprintf("%.1f", expansion_box[5].to_f)))
1363   }
1364   if h[:status] == 0
1365     @expansion_box[0] = h["xmin"].to_f
1366     @expansion_box[1] = h["xmax"].to_f
1367     @expansion_box[2] = h["ymin"].to_f
1368     @expansion_box[3] = h["ymax"].to_f
1369     @expansion_box[4] = h["zmin"].to_f
1370     @expansion_box[5] = h["zmax"].to_f
1371   else
1372     return
1373   end
1374
1375   #  Enumerate all atoms within the expansion box
1376   syms = self.symmetries
1377   h = Hash.new
1378   xmin = @expansion_box[0]
1379   xmax = @expansion_box[1]
1380   ymin = @expansion_box[2]
1381   ymax = @expansion_box[3]
1382   zmin = @expansion_box[4]
1383   zmax = @expansion_box[5]
1384   xmin_d = xmin.floor
1385   xmax_d = xmax.floor - 1
1386   ymin_d = ymin.floor
1387   ymax_d = ymax.floor - 1
1388   zmin_d = zmin.floor
1389   zmax_d = zmax.floor - 1
1390   syms.each_with_index { |sym, i|
1391     each_atom { |ap|
1392       fr = sym * ap.fract_r
1393           dx = fr.x.floor
1394           dy = fr.y.floor
1395           dz = fr.z.floor
1396           fr.x -= dx
1397           fr.y -= dy
1398           fr.z -= dz
1399           symopi = i * 1000000 + (50 - dx) * 10000 + (50 - dy) * 100 + 50 - dz
1400           (h[symopi] ||= []).push(ap.index)
1401           xmin_d.upto(xmax_d) { |xd|
1402             next if fr.x + xd < xmin || fr.x + xd > xmax
1403             ymin_d.upto(ymax_d) { |yd|
1404               next if fr.y + yd < ymin || fr.y + yd > ymax
1405                   zmin_d.upto(zmax_d) { |zd|
1406                     next if fr.z + zd < zmin || fr.z + zd > zmax
1407                         next if xd == 0 && yd == 0 && zd == 0
1408                     symopid = symopi + xd * 10000 + yd * 100 + zd
1409                         (h[symopid] ||= []).push(ap.index)
1410                   }
1411                 }
1412           }
1413         }
1414   }
1415   h.keys.sort.each { |key|
1416     sym = key / 1000000
1417         dx = key % 1000000 / 10000 - 50
1418         dy = key % 10000 / 100 - 50
1419         dz = key % 100 - 50
1420         next if sym == 0 && dx == 0 && dy == 0 && dz == 0
1421 #       puts "[#{sym},#{dx},#{dy},#{dz}]: #{h[key].inspect}"
1422         expand_by_symmetry(IntGroup[h[key]], sym, dx, dy, dz)
1423   }
1424 end
1425
1426 def cmd_show_ortep
1427   mol = self
1428   begin
1429     tmp = create_temp_dir("ortep", mol.name)
1430   rescue
1431     error_message_box($!.to_s)
1432         return
1433   end
1434   tepexe = "#{ResourcePath}/ortep3/ortep3"
1435   if $platform == "win"
1436     tepexe += ".exe"
1437   end
1438   tepdata = []
1439   tepbounds = [0, 0, 400, 400]
1440   descs = {
1441         "Atoms"=>[
1442           ["all", "Shaded", "black", ""],
1443           ["Li-Ar", "Principals", "black", ""],
1444           ["C", "Boundary", "black", ""],
1445           ["H", "Fixed", "black", "0.1"]
1446         ],
1447         "Bonds"=>[
1448           ["all", "all", "4 Shades", "black"],
1449           ["Li-Ar", "Li-Ar", "3 Shades", "black"],
1450           ["C", "Li-Ar", "2 Shades", "black"],
1451           ["H", "all", "Open", "black"]
1452         ]
1453   }
1454   mol.open_auxiliary_window("Show ORTEP", :resizable=>true, :has_close_box=>true) {
1455     tepview = nil  #  Forward declaration
1456         tab = "Atoms"
1457         columns = {
1458           "Atoms"=>[["atom list", 65], ["type", 80], ["color", 40], ["radius", 50]],
1459           "Bonds"=>[["list1", 65], ["list2", 65], ["type", 80], ["color", 40]]
1460         }
1461         atom_types = ["Boundary", "Principals", "Axes", "Shaded", "Fixed"]
1462         bond_types = ["Open", "1 Shade", "2 Shades", "3 Shades", "4 Shades"]
1463         colors = ["black", "red", "green", "blue", "cyan", "magenta", "yellow"]
1464         #  ORTEP attributes
1465         get_ortep_attr = lambda {
1466           attr = Hash.new
1467           make_atom_list = lambda { |s, ln|
1468             ss = s.scan(/[-.0-9A-Za-z]+/)
1469                 # puts ss.inspect
1470                 ss.inject(IntGroup[]) { |r, it|
1471                   if it == "all"
1472                     r.add(mol.all)
1473                         next r
1474                   elsif it =~ /-|(\.\.)/
1475                     s0 = Regexp.last_match.pre_match
1476                         s1 = Regexp.last_match.post_match
1477                   else
1478                     s0 = s1 = it
1479                   end
1480                   if s0 =~ /^[0-9]+$/
1481                     #  Atomic numbers
1482                         s0 = s0.to_i
1483                         if s1 !~ /^[0-9]+$/ || (s1 = s1.to_i) < s0
1484                           raise "Bad atom list specification: #{s} in line #{ln}"
1485                         end
1486                         r.add(s0..s1)
1487                   else
1488                     #  Atomic symbols
1489                         s0 = Parameter.builtin.elements.find_index { |p| p.name == s0.capitalize }
1490                         s1 = Parameter.builtin.elements.find_index { |p| p.name == s1.capitalize }
1491                         if s0 == nil || s1 == nil
1492                           raise "Bad element symbol specification: #{s} in line #{ln}"
1493                         end
1494                         (s0..s1).each { |an|
1495                           r.add(mol.atom_group { |ap| ap.atomic_number == an } )
1496                         }
1497                   end
1498                   r
1499                 }
1500       }
1501           #  Atoms
1502           at = []
1503           descs["Atoms"].each_with_index { |d, idx|
1504                 list = make_atom_list.call(d[0], idx)
1505                 type = atom_types.find_index { |it| it == d[1] } || 0
1506                 col = (colors.find_index { |it| it == d[2] } || 0) + 1
1507                 rad = d[3]
1508                 at.push([list, type, col, rad])
1509           }
1510           attr["atoms"] = at
1511           #  Bonds
1512           bd = []
1513           descs["Bonds"].each_with_index { |d, idx|
1514                 list1 = make_atom_list.call(d[0], idx)
1515                 list2 = make_atom_list.call(d[1], idx)
1516                 type = bond_types.find_index { |it| it == d[2] } || 0
1517                 col = (colors.find_index { |it| it == d[3] } || 0) + 1
1518                 bd.push([list1, list2, type, col])
1519           }
1520           attr["bonds"] = bd
1521           attr
1522         }
1523     on_update_ortep = lambda { |it|
1524           #  ORTEP attributes
1525           attr = get_ortep_attr.call
1526           #  Create ORTEP input in the temporary directory
1527           open(tmp + "/TEP.IN", "w") { |fp| mol.export_ortep(fp, attr) }
1528           #  Run ORTEP
1529           cwd = Dir.pwd
1530           Dir.chdir(tmp)
1531           if FileTest.exist?("TEP001.PRN")
1532             File.unlink("TEP001.PRN")
1533           end
1534           pid = call_subprocess([tepexe], "Running ORTEP", nil, "NUL", "NUL")
1535           if FileTest.exist?("TEP001.PRN")
1536             File.rename("TEP001.PRN", "TEP001.ps")
1537           end
1538           Dir.chdir(cwd)
1539           if pid != 0
1540             msg = "ORTEP execution in #{tmp} failed with status #{pid}."
1541         message_box(msg, "ORTEP Failed", :ok, :warning)
1542           elsif !FileTest.exist?(tmp + "/TEP001.ps")
1543         msg = "ORTEP execution in #{tmp} failed with unknown error."
1544         message_box(msg, "ORTEP Failed", :ok, :warning)
1545       else
1546             open(tmp + "/TEP001.ps", "r") { |fp|
1547                   tepdata.clear
1548                   tepbounds = [100000, 100000, -100000, -100000]
1549                   fp.each_line { |ln|
1550                     ln.chomp!
1551                     x, y, c = ln.split
1552                         if ln =~ /setrgbcolor/
1553                           tepdata.push(["color", x.to_f, y.to_f, c.to_f])
1554                           next
1555                         elsif ln =~ /setgray/
1556                           x = x.to_f
1557                           tepdata.push(["color", x, x, x])
1558                           next
1559                         end
1560                         x = x.to_i
1561                         y = y.to_i
1562                         if c == "m"
1563                           tepdata.push([x, y])
1564                         elsif c == "l"
1565                           tepdata[-1].push(x, y)
1566                         else
1567                           next
1568                         end
1569                         tepbounds[0] = x if x < tepbounds[0]
1570                         tepbounds[1] = y if y < tepbounds[1]
1571                         tepbounds[2] = x if x > tepbounds[2]
1572                         tepbounds[3] = y if y > tepbounds[3]
1573                   }
1574                   #  [x, y, width, height]
1575                   tepbounds[2] -= tepbounds[0]
1576                   tepbounds[3] -= tepbounds[1]
1577                 }
1578                 tepview.refresh_rect(tepview[:frame])
1579           end
1580         }
1581         on_draw_ortep = lambda { |it|
1582           clear
1583           frame = it[:frame]
1584           brush(:color=>[1, 1, 1, 1])
1585           pen(:color=>[1, 1, 1, 1])
1586           draw_rectangle(frame)
1587           pen(:color=>[0, 0, 0, 1])
1588           rx = (frame[2] - 10.0) / tepbounds[2]
1589           ry = (frame[3] - 10.0) / tepbounds[3]
1590           rx = ry if rx > ry
1591           dx = (frame[2] - tepbounds[2] * rx) * 0.5
1592           dy = (frame[3] + tepbounds[3] * rx) * 0.5
1593           tepdata.each { |d|
1594             if d[0] == "color"
1595                   pen(:color=>[d[1], d[2], d[3], 1])
1596                   next
1597                 end
1598             x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
1599                 y0 = (dy - (d[1] - tepbounds[1]) * rx).to_i
1600                 (d.length / 2 - 1).times { |i|
1601                   x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
1602                   y1 = (dy - (d[i * 2 + 3] - tepbounds[1]) * rx).to_i
1603                   draw_line(x0, y0, x1, y1)
1604                   x0 = x1
1605                   y0 = y1
1606                 }
1607           }
1608         }
1609         on_export_eps = lambda { |fname|
1610           frame = item_with_tag("ortep")[:frame].dup
1611           dx = dy = 5
1612           rx = (frame[2] - dx * 2) / tepbounds[2]
1613           ry = (frame[3] - dy * 2) / tepbounds[3]
1614           rx = ry if rx > ry
1615           frame[2] = (rx * tepbounds[2] + dx * 2).to_i
1616           frame[3] = (rx * tepbounds[3] + dy * 2).to_i
1617           #  Assumes A4 paper and center the bounding box
1618           frame[0] = (595 - frame[2]) / 2
1619           frame[1] = (842 - frame[3]) / 2
1620           xmax = frame[0] + frame[2]
1621           ymax = frame[1] + frame[3]
1622           open(fname, "w") { |fp|
1623             fp.print "%!PS-Adobe-3.0 EPSF-3.0
1624 %%Creator: Molby
1625 %%BoundingBox: #{frame[0]} #{frame[1]} #{xmax} #{ymax}
1626 %%Pages: 1
1627 %%Orientation: Landscape
1628 %%BeginProlog
1629 /m {moveto} def
1630 /l {lineto} def
1631 %%EndProlog
1632 %%Page: 1 1
1633 %%BeginPageSetup
1634 0 setgray 1 setlinecap 1 setlinewidth
1635 %%EndPageSetup
1636 "
1637             tepdata.each { |d|
1638               if d[0] == "color"
1639                     fp.printf "%.3f %.3f %.3f setrgbcolor\n", d[1], d[2], d[3]
1640                     next
1641                   end
1642               x0 = frame[0] + dx + (d[0] - tepbounds[0]) * rx
1643                   y0 = frame[1] + dy + (d[1] - tepbounds[1]) * rx
1644                   fp.printf "%.2f %.2f m\n", x0, y0
1645                   (d.length / 2 - 1).times { |i|
1646                     x1 = frame[0] + dx + (d[i * 2 + 2] - tepbounds[0]) * rx
1647                     y1 = frame[1] + dy + (d[i * 2 + 3] - tepbounds[1]) * rx
1648                     fp.printf "%.2f %.2f l\n", x1, y1
1649                     x0 = x1
1650                     y0 = y1
1651                   }
1652                   fp.print "stroke\n"
1653             }
1654         fp.print "showpage\n%%Trailer\n%%EOF\n"
1655           }
1656         }
1657         on_export_bitmap = lambda { |fname|
1658           frame = item_with_tag("ortep")[:frame].dup
1659           dx = dy = 5
1660           scale = 5.0
1661           rx = (frame[2] - dx * 2) * scale / tepbounds[2]
1662           ry = (frame[3] - dy * 2) * scale / tepbounds[3]
1663           rx = ry if rx > ry
1664           frame[2] = (rx * tepbounds[2] + dx * 2).to_i
1665           frame[3] = (rx * tepbounds[3] + dy * 2).to_i
1666       bmp = Bitmap.new(frame[2], frame[3], 32)
1667           bmp.focus_exec {
1668             clear
1669             pen(:color=>[0, 0, 0, 1], :width=>scale)
1670             tepdata.each { |d|
1671               if d[0] == "color"
1672                     pen(:color=>[d[1], d[2], d[3], 1])
1673                     next
1674                   end
1675               x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
1676                   y0 = (dy + (tepbounds[3] - d[1] + tepbounds[1]) * rx).to_i
1677                   (d.length / 2 - 1).times { |i|
1678                     x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
1679                     y1 = (dy + (tepbounds[3] - d[i * 2 + 3] + tepbounds[1]) * rx).to_i
1680                     draw_line(x0, y0, x1, y1)
1681                     x0 = x1
1682                     y0 = y1
1683                   }
1684             }
1685           }
1686           bmp.save_to_file(fname)
1687         }
1688         on_export = lambda { |it|
1689           basename = (mol.path ? File.basename(mol.path, ".*") : mol.name)
1690       fname = Dialog.save_panel("Export ORTEP file:", mol.dir, basename + ".eps",
1691             "Encapsulated PostScript (*.eps)|*.eps|PNG (*.png)|*.png|TIFF (*.tif)|*.tif|ORTEP Instruction (*.tep)|*.tep|All files|*.*")
1692           return if !fname
1693           ext = File.extname(fname).downcase
1694           if ext == ".eps"
1695             on_export_eps.call(fname)
1696           elsif ext == ".png" || ext == ".tif" || ext == ".tiff"
1697             on_export_bitmap.call(fname)
1698           elsif ext == ".tep"
1699             filecopy("#{tmp}/TEP.IN", fname)
1700       end
1701         }
1702         @on_document_modified = lambda { |m|
1703         }
1704         #  Close handler (called when the close box is pressed or the document is closed)
1705         @on_close = lambda { |*d|
1706           erase_old_logs(tmp)
1707           tmp = nil
1708           true
1709         }
1710         on_select_tab = lambda { |it|
1711           tab = ["Atoms", "Bonds"][it[:value]]
1712           # puts "tab = #{tab}"
1713           table = item_with_tag("table")
1714           table[:columns] = columns[tab]
1715           table[:selection] = IntGroup[]
1716           table[:refresh] = true
1717         }
1718         get_count = lambda { |it| descs[tab].count }
1719         get_value = lambda { |it, row, col| descs[tab][row][col] }
1720         has_popup_menu = lambda { |it, row, col|
1721           if tab == "Atoms"
1722             if col == 1
1723               return atom_types
1724                 elsif col == 2
1725                   return colors
1726                 end
1727           elsif tab == "Bonds"
1728                 if col == 2
1729               return bond_types
1730                 elsif col == 3
1731                   return colors
1732                 end
1733           end
1734           return nil
1735         }
1736         on_selection_changed = lambda { |it|
1737           sel = it[:selection]
1738           item_with_tag("remove")[:enabled] = (sel == nil || sel.count == 0 ? false : true)
1739         }
1740         is_item_editable = lambda { |it, row, col|
1741           return has_popup_menu.call(it, row, col) == nil
1742         }
1743         popup_menu_selected = lambda { |it, row, col, sel|
1744           titles = has_popup_menu.call(it, row, col)
1745           return nil if titles == nil
1746           descs[tab][row][col] = titles[sel]
1747           it[:refresh] = true
1748         }
1749         set_value = lambda { |it, row, col, val|
1750           descs[tab][row][col] = val
1751           it[:refresh] = true
1752         }
1753         add_to_table = lambda { |it|
1754           table = item_with_tag("table")
1755           d = descs[tab][-1].dup
1756           descs[tab].push(d)
1757           table[:refresh] = true
1758           table[:selection] = IntGroup[descs[tab].count - 1]
1759           # puts descs[tab].inspect
1760         }
1761         remove_from_table = lambda { |it|
1762           table = item_with_tag("table")
1763           sel = table[:selection]
1764           sel.reverse_each { |n|
1765             descs[tab][n, 1] = []
1766           }
1767           table[:refresh] = true
1768         }
1769         layout(2,
1770           layout(1,
1771             item(:popup, :subitems=>["Atoms", "Bonds"], :action=>on_select_tab),
1772                 item(:table, :width=>240, :height=>380, :flex=>[0,0,0,0,1,1], :tag=>"table",
1773                   :columns=>columns["Atoms"],
1774                   :on_count=>get_count,
1775                   :on_get_value=>get_value,
1776                   :on_set_value=>set_value,
1777                   :on_selection_changed=>on_selection_changed,
1778                   :is_item_editable=>is_item_editable,
1779                   :has_popup_menu=>has_popup_menu,
1780                   :on_popup_menu_selected=>popup_menu_selected),
1781                 layout(2,
1782                   item(:button, :title=>"Add", :tag=>"add", :action=>add_to_table),
1783                   item(:button, :title=>"Remove", :tag=>"remove", :action=>remove_from_table, :enabled=>false),
1784                   :flex=>[0,1,0,0,0,0]),
1785             :flex=>[0,0,0,0,0,1]),
1786           layout(1,
1787             item(:view, :width=>400, :height=>400, :tag=>"ortep", :on_paint=>on_draw_ortep, :flex=>[0,0,0,0,1,1]),
1788             layout(2,
1789               item(:button, :title=>"Refresh", :action=>on_update_ortep, :align=>:left),
1790                   item(:button, :title=>"Export as...", :action=>on_export, :align=>:right),
1791                   :flex=>[0,1,0,0,0,0]),
1792             :flex=>[0,0,0,0,1,1]),
1793           :flex=>[0,0,0,0,1,1]
1794         )
1795         set_min_size(480, 200)
1796         tepview = item_with_tag("ortep")
1797 #       listen(mol, "documentModified", on_document_modified)
1798 #       listen(mol, "documentWillClose", lambda { |m| close } )
1799         @on_document_modified.call(nil)
1800     on_update_ortep.call(nil)
1801         show
1802   }
1803 end
1804
1805   def cmd_unit_cell
1806     mol = self
1807         new_box = nil
1808     hash = Dialog.run("Unit Cell") {
1809           @mol = mol
1810           @box = @mol.box
1811           @box = (@box ? @box.dup : nil)
1812           @box_save = nil
1813
1814           def set_cell_value(item1)
1815                 alpha = value("alpha").to_f * PI / 180.0
1816                 beta = value("beta").to_f * PI / 180.0
1817                 gamma = value("gamma").to_f * PI / 180.0
1818                 if alpha.abs < 1e-8 || beta.abs < 1e-8 || gamma.abs < 1e-8
1819                   return
1820                 end
1821             if @box == nil
1822                   @box = [Vector3D[1,0,0],Vector3D[0,1,0],Vector3D[0,0,1],Vector3D[0,0,0],[1,1,1]]
1823                 end
1824                 av = Vector3D[1, 0, 0]
1825                 bv = Vector3D[cos(gamma), sin(gamma), 0]
1826                 cx = cos(beta)
1827                 cy = (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma)
1828                 cz = sqrt_safe(1.0 - cx * cx - cy * cy)
1829                 cv = Vector3D[cx, cy, cz]
1830                 x0 = @box[0].normalize
1831                 z0 = (@box[0].cross(@box[1])).normalize
1832                 y0 = (z0.cross(x0)).normalize
1833                 tr = Transform[x0, y0, z0, Vector3D[0, 0, 0]]
1834                 av = (tr * av) * value("a").to_f
1835                 bv = (tr * bv) * value("b").to_f
1836                 cv = (tr * cv) * value("c").to_f
1837                 @box[0] = av
1838                 @box[1] = bv
1839                 @box[2] = cv
1840                 3.times { |i|
1841                   3.times { |j|
1842                     if @box[i][j].abs < 1e-8
1843                           @box[i][j] = 0.0
1844                         end
1845                   }
1846                 }
1847                 update_values(item1)
1848           end
1849
1850           def set_box_value(item1)
1851             h = Hash.new
1852                 enabled = false
1853                 ["o0", "o1", "o2", "a0", "a1", "a2", "b0", "b1", "b2", "c0", "c1", "c2"].each { |k|
1854                   begin
1855                     s = value(k)
1856                         if s == nil || s == ""
1857                           h[k] = 0.0
1858                         else
1859                           enabled = true
1860                       h[k] = Float(eval(s))
1861                           set_value(k, h[k].to_s)
1862                         end
1863                   rescue
1864                     mes = "Cannot evaluate #{value(k)}: " + $!.to_s
1865                         Dialog.run("Value Error") {
1866                           layout(1, item(:text, :title=>mes))
1867                           set_attr(1, :hidden=>true)
1868                         }
1869                         return nil
1870                   end
1871                 }
1872                 if enabled
1873                   ax = Vector3D[h["a0"], h["a1"], h["a2"]]
1874                   bx = Vector3D[h["b0"], h["b1"], h["b2"]]
1875                   cx = Vector3D[h["c0"], h["c1"], h["c2"]]
1876                   ox = Vector3D[h["o0"], h["o1"], h["o2"]]
1877                   fx = [1, 1, 1]
1878                   if ax.length2 < 1e-8
1879                     fx[0] = 0
1880                     ax = Vector3D[1, 0, 0]
1881                   end
1882                   if bx.length2 < 1e-8
1883                     fx[1] = 0
1884                     bx = Vector3D[0, 1, 0]
1885                   end
1886                   if cx.length2 < 1e-8
1887                     fx[2] = 0
1888                     cx = Vector3D[0, 0, 1]
1889                   end
1890                   @box = [ax, bx, cx, ox, fx]
1891                 else
1892                   @box = nil
1893                 end
1894                 update_values(item1)
1895           end
1896           def clear_box(item1)
1897             @box_save = @box
1898                 @box = nil
1899                 set_attr("restore", :enabled=>true)
1900                 update_values(item1)
1901           end
1902           def restore_box(item1)
1903             @box = @box_save.dup
1904                 set_attr("restore", :enabled=>false)
1905                 update_values(item1)
1906           end
1907           def angle(v1, v2)
1908             acos(v1.dot(v2)/(v1.length * v2.length)) * 180.0 / PI
1909           end
1910           def update_values(it)
1911             ["a0", "a1", "a2", "b0", "b1", "b2", "c0", "c1", "c2", "o0", "o1", "o2",
1912                  "a", "b", "c", "alpha", "beta", "gamma"].each_with_index { |tag, i|
1913                   if @box == nil
1914                     set_value(tag, "")
1915                   else
1916                     if i < 12
1917                       val = @box[i / 3][i % 3]
1918                           fmt = "%.8g"
1919                     elsif i < 15
1920                       val = @box[i - 12].length
1921                           fmt = "%g"
1922                     else
1923                       val = angle(@box[(i - 14) % 3], @box[(i - 13) % 3])
1924                           fmt = "%g"
1925                         end
1926                         set_value(tag, sprintf(fmt, val))
1927                   end
1928                 }
1929           end
1930
1931           layout(4,
1932             item(:text, :title=>"Cell parameters:"),
1933                 -1, -1, -1,
1934
1935                 item(:text, :title=>"a/b/c"),
1936                 item(:textfield, :width=>140, :tag=>"a"),
1937                 item(:textfield, :width=>140, :tag=>"b"),
1938                 item(:textfield, :width=>140, :tag=>"c"),
1939
1940                 item(:text, :title=>"α/β/γ"),
1941                 item(:textfield, :width=>140, :tag=>"alpha"),
1942                 item(:textfield, :width=>140, :tag=>"beta"),
1943                 item(:textfield, :width=>140, :tag=>"gamma"),
1944
1945                 -1, -1, -1,
1946         item(:button, :title=>"Set Cell Parameters", :align=>:right, :action=>:set_cell_value),
1947         item(:line),
1948                 -1, -1, -1,
1949
1950             item(:text, :title=>"Axis vectors:"),
1951                 -1, -1, -1,
1952
1953             item(:text, :title=>"origin"),
1954                 item(:textfield, :width=>140, :tag=>"o0"),
1955                 item(:textfield, :width=>140, :tag=>"o1"),
1956                 item(:textfield, :width=>140, :tag=>"o2"),
1957
1958             item(:text, :title=>"a-axis"),
1959                 item(:textfield, :width=>140, :tag=>"a0"),
1960                 item(:textfield, :width=>140, :tag=>"a1"),
1961                 item(:textfield, :width=>140, :tag=>"a2"),
1962
1963             item(:text, :title=>"b-axis"),
1964                 item(:textfield, :width=>140, :tag=>"b0"),
1965                 item(:textfield, :width=>140, :tag=>"b1"),
1966                 item(:textfield, :width=>140, :tag=>"b2"),
1967
1968             item(:text, :title=>"c-axis"),
1969                 item(:textfield, :width=>140, :tag=>"c0"),
1970                 item(:textfield, :width=>140, :tag=>"c1"),
1971                 item(:textfield, :width=>140, :tag=>"c2"),
1972
1973         -1, -1, -1,
1974         item(:button, :title=>"Set Axis Vectors", :align=>:right, :action=>:set_box_value),
1975         item(:line),
1976                 -1, -1, -1,
1977
1978                 item(:text, :title=>"(Ruby expressions are allowed as the values; e.g. 12.0*cos(60*PI/180))"),
1979                 -1, -1, -1,
1980
1981                 item(:button, :title=>"Clear Cell", :tag=>"clear", :action=>:clear_box),
1982                 item(:button, :title=>"Restore Cell", :tag=>"restore", :action=>:restore_box, :enabled=>false),
1983                 -1, -1)
1984           set_attr(0, :action=>lambda { |it|
1985             if set_box_value(it)
1986                   new_box = @box
1987                   end_modal(it)
1988                 end
1989           } )
1990           update_values(nil)
1991         }
1992         if hash[:status] == 0
1993           if new_box != nil
1994             h = Dialog.run(" ") {
1995                   layout(1,
1996                     item(:text, :title=>"Do you want to keep the fractional coordinates?"),
1997                         item(:radio, :title=>"Yes, convert the cartesian coordinates so that the fractional coordinates will be unchanged.",
1998                           :tag=>"yes", :value=>1),
1999                         item(:radio, :title=>"No, keep the cartesian coordinates unchanged.", :tag=>"no"))
2000                   radio_group("yes", "no")
2001                 }
2002                 if h[:status] == 0
2003                   yes = h["yes"] != 0
2004                   new_box.push(yes)
2005                   mol.box = new_box
2006                 end
2007           else
2008             mol.box = nil
2009           end
2010         end
2011   end
2012
2013
2014 #  For debug; check the symmetry data in space_groups.txt for self-consistency
2015 def Molecule.check_space_group
2016   zero_transform = Transform.zero
2017   group_member_p = lambda { |g, tr|
2018         g.each_with_index { |tr2, idx|
2019           tr2 = tr2 - tr
2020           #  Wrap the translational part to [-0.5..0.5]
2021           3.times { |k| tr2[3, k] -= (tr2[3, k] + 0.5).floor }
2022           if tr2.nearly_equal?(zero_transform)
2023                 return idx
2024           end
2025         }
2026         return nil
2027   }
2028   @@space_groups.each { |g|
2029     name = g[0]
2030         group = g[1]
2031         n = group.count
2032         err = ""
2033         n.times { |i|
2034           tri = group[i]
2035           begin
2036             trx = tri.inverse
2037           rescue
2038             err += "The element #{i + 1} is not regular transform; #{trx ? symmetry_to_string(trx) : 'nil'}\n"
2039                 next
2040           end
2041           3.times { |k| trx[3, k] -= trx[3, k].floor }  #  Wrap the translation part to [0..1]
2042           if !group_member_p.call(group, trx)
2043                 err += "The element #{i + 1} has no inverse element; #{symmetry_to_string(trx)}\n"
2044                 next
2045           end
2046           n.times { |j|
2047                 trj = group[j]
2048                 trx = tri * trj
2049                 3.times { |k| trx[3, k] -= trx[3, k].floor }  #  Wrap the translation part to [0..1]
2050                 if !group_member_p.call(group, trx)
2051                   err += "The element #{i + 1} * #{j + 1} is not in the group; "
2052                   err += "#{symmetry_to_string(tri)} * #{symmetry_to_string(trj)} = #{symmetry_to_string(trx)}\n"
2053                 end
2054           }
2055         }
2056         if err == ""
2057           msg = "#{name} is a complete group\n"
2058         else
2059           msg = "#{name} is NOT a complete group\n" + err
2060           error_message_box(msg)
2061         end
2062 #       break if err != ""
2063   }
2064 end
2065
2066 def Molecule.setup_space_groups
2067   if defined?(@@space_groups)
2068     return @@space_groups
2069   end
2070
2071   #  The debug flag causes self-check of the space group table
2072   debug = false
2073
2074   lines = [] if debug
2075
2076   @@space_groups = []       #  Sequential table of space groups (array of [name, array_of_transform])
2077   @@space_group_list = []   #  [["Triclinic",    -> Crystal system
2078                                                         #     ["#1: P1",     -> Space group family
2079                                                         #       [["P1", 0]]],-> Variations
2080                                                         #     ...],
2081                                                         #   ["Monoclinic",
2082                                                         #     ...
2083                                                         #     ["#7: Pc",
2084                                                         #       [["Pc", 6], ["Pa", 7], ["Pn", 8]]],
2085                                                         #     ...]]
2086   last_group_code = 0
2087   Kernel.open(ScriptPath + "/space_groups.txt", "r") { |fp|
2088         while ln = fp.gets
2089           lines.push(ln) if debug
2090           if ln =~ /\[(\w+)\]/
2091                 label = $1
2092                 @@space_group_list.push([label, []])
2093                 next
2094           end
2095           if ln =~ /^(\d+) +(.*)$/
2096                 group_code = $1.to_i
2097                 group_name = $2
2098             name = "#" + $1 + ": " + $2
2099                 group = []
2100                 group_gen = []
2101                 while ln = fp.gets
2102                   lines.push(ln.dup) if debug
2103                   ln.chomp!
2104                   if ln =~ /^\s*$/
2105                         break
2106                   elsif ln =~ /^\+\(/
2107                         lattice = ln.scan(/((^\+)|,)\(([\/,0-9]+)\)/).map {
2108                           |a| a[2].split(/, */).map {
2109                                 |x| if x =~ /^(\d+)\/(\d+)$/
2110                                       $1.to_f / $2.to_f
2111                                         else
2112                                           x.to_f
2113                                         end } }
2114                         lattice.each { |lat|
2115                           group.each { |tr|
2116                                 tr = tr.dup
2117                                 tr[3, 0] += lat[0]
2118                                 tr[3, 1] += lat[1]
2119                                 tr[3, 2] += lat[2]
2120                                 group_gen.push(tr)
2121                           }
2122                         }
2123                   else
2124                         #  Rectify the shorthand descriptions
2125                         ln1 = ln.dup if debug
2126                         mod = false
2127                         mod = ln.gsub!(/(\d)(\d)/, "\\1/\\2")
2128                         mod = ln.gsub!(/(^|[^\/])(\d)([^\/]|$)/, "\\11/\\2\\3") || mod
2129                         mod = ln.gsub!(/([0-9xyzXYZ])([xyzXYZ])/, "\\1+\\2") || mod
2130                         if debug && mod
2131                           puts "#{fp.lineno}: #{ln1} -> #{ln}\n"
2132                         end
2133                         begin
2134                           group.push(string_to_symmetry(ln))
2135                           if debug
2136                                 lines[-1] = symmetry_to_string(group[-1]).gsub(/ +/, "") + "\n"
2137                           end
2138                         rescue
2139                           error_message_box($!.message + "(line #{fp.lineno})")
2140                         end
2141                   end
2142                 end
2143                 group.concat(group_gen)
2144                 @@space_groups.push([name, group])
2145                 group_index = @@space_groups.count - 1
2146                 if last_group_code != group_code
2147                   idx = (name.index(" (") || 0) - 1  #  Remove parenthesized phrase
2148                   group_family_name = name[0..idx]
2149                   @@space_group_list[-1][1].push([group_family_name, []])
2150                   last_group_code = group_code
2151                 end
2152                 @@space_group_list[-1][1][-1][1].push([group_name, group_index])
2153           end
2154         end
2155   }
2156   if debug
2157         Molecule.check_space_group
2158         fn = Dialog.save_panel("Save modified space_group.txt:")
2159         if fn
2160           Kernel.open(fn, "w") { |fp|
2161                 lines.each { |ln|
2162                   fp.print ln
2163                 }
2164           }
2165         end
2166   end
2167   return @@space_groups
2168 end
2169
2170 def cmd_symmetry_operations
2171
2172   mol = self
2173   syms = mol.symmetries
2174
2175   if !defined?(@@space_groups)
2176     Molecule.setup_space_groups
2177   end
2178
2179   #  Menu titles
2180   crystal_systems = @@space_group_list.map { |m| m[0] }
2181   space_groups = @@space_group_list[0][1].map { |m| m[0] }
2182   variations = @@space_group_list[0][1][0][1].map { |m| m[0] }
2183
2184   #  The index (for @@space_groups) of the current space group (an integer or nil)
2185   current_space_group = nil
2186
2187   #  Find current space group
2188   guess_space_group = lambda {
2189     current_space_group = nil
2190     s = syms.map { |tr| tr = tr.dup; 3.times { |k| tr[3, k] -= tr[3, k].floor }; tr }
2191         @@space_groups.each_with_index { |g, i|
2192           next if g[1].count != s.count
2193           ss = s.dup
2194           g[1].each { |tr|
2195             idx = ss.find_index { |tr2| tr.nearly_equal?(tr2) }
2196                 break if idx == nil
2197                 ss.delete_at(idx)
2198           }
2199           if ss.empty?
2200             current_space_group = i
2201                 break
2202           end
2203         }
2204   }
2205
2206   select_space_group = lambda { |it|
2207
2208     h = Dialog.run("Select Space Group") {
2209           update_operation = lambda { |*d|
2210             it = item_with_tag("operations")
2211                 cval = item_with_tag("crystal_system")[:value]
2212                 sval = item_with_tag("space_group")[:value]
2213                 vval = item_with_tag("variation")[:value]
2214                 idx = @@space_group_list[cval][1][sval][1][vval][1]
2215             op = ""
2216             @@space_groups[idx][1].each { |tr|
2217               op += symmetry_to_string(tr) + "\n"
2218             }
2219             it[:value] = op
2220           }
2221       crystal_system_popup_handler = lambda { |it|
2222             val = it[:value]
2223         space_groups = @@space_group_list[val][1].map { |m| m[0] }
2224         variations = @@space_group_list[val][1][0][1].map { |m| m[0] }
2225                 item_with_tag("space_group")[:subitems] = space_groups
2226                 item_with_tag("space_group")[:value] = 0
2227                 item_with_tag("variation")[:subitems] = variations
2228                 item_with_tag("variation")[:value] = 0
2229                 update_operation.call
2230       }
2231       space_group_popup_handler = lambda { |it|
2232             val = it[:value]
2233                 cval = item_with_tag("crystal_system")[:value]
2234                 variations = @@space_group_list[cval][1][val][1].map { |m| m[0] }
2235                 item_with_tag("variation")[:subitems] = variations
2236                 item_with_tag("variation")[:value] = 0
2237                 update_operation.call
2238       }
2239           variation_popup_handler = lambda { |it|
2240             update_operation.call
2241           }
2242           layout(2,
2243             item(:text, :title=>"Crystal System"),
2244                 item(:popup, :subitems=>crystal_systems, :tag=>"crystal_system", :action=>crystal_system_popup_handler),
2245                 item(:text, :title=>"Space Group"),
2246                 item(:popup, :subitems=>space_groups, :tag=>"space_group", :action=>space_group_popup_handler),
2247                 item(:text, :title=>"Variation"),
2248                 item(:popup, :subitems=>variations, :tag=>"variation", :action=>variation_popup_handler),
2249                 item(:textview, :width=>360, :height=>200, :editable=>false, :tag=>"operations"),
2250                 -1)
2251           if current_space_group
2252             catch(:end_loop) {
2253               @@space_group_list.each_with_index { |m, i|
2254                     m[1].each_with_index { |mm, j|
2255                       mm[1].each_with_index { |mmm, k|
2256                             if mmm[1] == current_space_group
2257                               (it = item_with_tag("crystal_system"))[:value] = i
2258                                   crystal_system_popup_handler.call(it)
2259                               (it = item_with_tag("space_group"))[:value] = j
2260                                   space_group_popup_handler.call(it)
2261                               (it = item_with_tag("variation"))[:value] = k
2262                                   variation_popup_handler.call(it)
2263                                   throw(:end_loop)
2264                             end
2265                           }
2266                     }
2267                   }
2268                 }
2269           end
2270         }
2271         if h[:status] == 0
2272           cval = h["crystal_system"]
2273           sval = h["space_group"]
2274           vval = h["variation"]
2275           idx = @@space_group_list[cval][1][sval][1][vval][1]
2276           syms = @@space_groups[idx][1].dup
2277         end
2278   }
2279
2280   h = Dialog.run("Symmetry Operations") {
2281     update_space_group_text = lambda {
2282           guess_space_group.call
2283           if current_space_group
2284             label = @@space_groups[current_space_group][0]
2285           else
2286             label = ""
2287           end
2288           item_with_tag("space_group")[:value] = label
2289         }
2290     layout(1,
2291           layout(2,
2292             item(:text, :title=>"Space Group:"),
2293                 item(:textfield, :width=>200, :editable=>false, :tag=>"space_group")),
2294           item(:button, :title=>"Select...", :action=>lambda { |it|
2295                 select_space_group.call(it)
2296                 item_with_tag("sym_table")[:refresh] = true
2297                 update_space_group_text.call }, :align=>:right),
2298           item(:table, :tag=>"sym_table", :width=>300, :height=>300,
2299                 :columns=>[["Symmetry Operations", 280]],
2300         :on_count=>lambda { |it| syms.count },
2301         :on_get_value=>lambda { |it, row, col| symmetry_to_string(syms[row]) },
2302         :on_set_value=>lambda { |it, row, col, val|
2303                   syms[row] = string_to_symmetry(val)
2304                   update_space_group_text.call },
2305         :is_item_editable=>lambda { |it, row, col| true },
2306         :on_selection_changed=>lambda { |it|
2307                   item_with_tag("delete")[:enabled] = (it[:selection].count > 0) } ),
2308           layout(2,
2309             item(:togglebutton, :title=>"+", :width=>40, :tag=>"add",
2310                   :action=>lambda { |it|
2311                     syms.push(Transform.new); item_with_tag("sym_table")[:refresh] = true
2312                         update_space_group_text.call } ),
2313                 item(:togglebutton, :title=>"-", :width=>40, :tag=>"delete", :enabled=>false,
2314                   :action=>lambda { |it|
2315                     item_with_tag("sym_table")[:selection].reverse_each { |n| syms.delete_at(n) }
2316                         item_with_tag("sym_table")[:refresh] = true
2317                         update_space_group_text.call } ),
2318                 :padding=>0, :margin=>0 ) )
2319         item_with_tag("sym_table")[:refresh] = true
2320         update_space_group_text.call
2321   }
2322   if h[:status] == 0
2323     self.remove_symmetries(nil)
2324         syms.each { |s|
2325           self.add_symmetry(s)
2326     }
2327   end
2328 end
2329
2330 def cmd_show_periodic_image
2331     mol = self
2332     hash = Dialog.run("Show Periodic Image") {
2333           @mol = mol
2334           def set_periodic_image(it)
2335             a = []
2336             ["amin", "amax", "bmin", "bmax", "cmin", "cmax"].each_with_index { |k, i|
2337                   s = value(k)
2338                   if s == nil || s == ""
2339                     a[i] = 0
2340                   else
2341                     a[i] = Integer(s)
2342                   end
2343                 }
2344             @mol.show_periodic_image(a)
2345                 @mol.show_periodic_image = (attr("show_flag", :value) == 1 ? true : false)
2346           end
2347           pimage = @mol.show_periodic_image
2348           flag = @mol.show_periodic_image?
2349           layout(4,
2350             item(:checkbox, :title=>"Show Periodic Image", :tag=>"show_flag", :value=>(flag ? 1 : 0),
2351                   :action=>lambda { |it| @mol.show_periodic_image = (it[:value] == 1 ? true : false) } ),
2352                 -1, -1, -1,
2353                 item(:text, :title=>"a-axis"),
2354                 item(:textfield, :width=>80, :tag=>"amin", :value=>pimage[0].to_s),
2355                 item(:text, :title=>"to"),
2356                 item(:textfield, :width=>80, :tag=>"amax", :value=>pimage[1].to_s),
2357                 item(:text, :title=>"b-axis"),
2358                 item(:textfield, :width=>80, :tag=>"bmin", :value=>pimage[2].to_s),
2359                 item(:text, :title=>"to"),
2360                 item(:textfield, :width=>80, :tag=>"bmax", :value=>pimage[3].to_s),
2361                 item(:text, :title=>"c-axis"),
2362                 item(:textfield, :width=>80, :tag=>"cmin", :value=>pimage[4].to_s),
2363                 item(:text, :title=>"to"),
2364                 item(:textfield, :width=>80, :tag=>"cmax", :value=>pimage[5].to_s),
2365                 item(:button, :title=>"Set", :action=>:set_periodic_image))
2366           set_attr(0, :action=>lambda { |it| set_periodic_image(it); end_modal(it) } )
2367         }
2368 end
2369
2370 def has_expanded_atoms
2371   atoms.find { |ap| ap.symop != nil }
2372 end
2373
2374 def cmd_remove_expanded_atoms
2375   mol = self
2376   return unless mol.has_expanded_atoms
2377   hash = Dialog.run("Remove Expanded Atoms") {
2378     layout(1,
2379           item(:radio, :title=>"Remove all expanded atoms.", :tag=>"all", :value=>1),
2380           item(:radio, :title=>"Keep the expanded atoms that are in the same fragment with original atoms.", :tag=>"partial"))
2381         radio_group("all", "partial")
2382   }
2383   if hash[:status] == 0
2384     if hash["all"] == 1
2385           rm = mol.atom_group { |ap| ap.symop != nil }
2386         else
2387           rm = IntGroup[]
2388           mol.each_fragment { |f|
2389             if !f.find { |i| mol.atoms[i].symop == nil }
2390                   rm.add(f)
2391                 end
2392           }
2393         end
2394         mol.remove(rm)
2395   end
2396 end
2397
2398 end
2399
2400 require_cell = lambda { |m| m && m.cell != nil }
2401
2402 register_menu("Xtal\tUnit Cell...", :cmd_unit_cell)
2403 register_menu("Xtal\tSymmetry Operations...", :cmd_symmetry_operations)
2404 register_menu("Xtal\t-", nil)
2405 register_menu("Xtal\tComplete by Symmetry", :complete_by_symmetry, lambda { |m| m && (m.cell != nil || m.nsymmetries > 1) } )
2406 register_menu("Xtal\tCreate Packing Diagram...", :create_packing_diagram, require_cell)
2407 register_menu("Xtal\tShow Periodic Image...", :cmd_show_periodic_image, require_cell)
2408 register_menu("Xtal\tRemove Expanded Atoms...", :cmd_remove_expanded_atoms, lambda { |m| m && m.has_expanded_atoms } )
2409 register_menu("Xtal\t-", nil)
2410 register_menu("Xtal\tBonds and Angles with Sigma...", :cmd_bond_angle_with_sigma, :non_empty)
2411 register_menu("Xtal\tBest-fit Planes...", :cmd_plane, :non_empty)
2412 register_menu("Xtal\tShow ORTEP...", :cmd_show_ortep, :non_empty)