# coding: utf-8
# 
#  crystal.rb
#
#  Created by Toshi Nagata.
#  Copyright 2012 Toshi Nagata. All rights reserved.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation version 2 of the License.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

#  Definition for use in ORTEP
class AtomRef
  def to_adc
    sym = self.symop
    if sym == nil
      idx = self.index + 1
      symcode = 55501
    else
      idx = sym[4] + 1
      symcode = (sym[1] + 5) * 10000 + (sym[2] + 5) * 1000 + (sym[3] + 5) * 100 + sym[0] + 1
    end
    return idx, symcode
  end
end

module Kernel

def symmetry_to_string(sym)
  #  Sym is a Transform object defined in crystallographic coordinates
  str = ""
  3.times { |j|
    s = ""
	a = sym[3, j]
	if a.abs > 1e-4
	  if a > 0.0 && s != ""
	    s += "+"
	  elsif a < 0.0
	    s += "-"
	  end
	  a = a.abs
	  if (a - (n = (a + 0.5).floor)).abs < 1e-4
	    #  integer
		s += n.to_s
	  elsif ((a * 2) - (n = (a * 2 + 0.5).floor)).abs < 1e-4
	    #  n/2
	    s += n.to_s + "/2"
	  elsif ((a * 3) - (n = (a * 3 + 0.5).floor)).abs < 1e-4
	    #  n/3
		s += n.to_s + "/3"
	  elsif ((a * 4) - (n = (a * 4 + 0.5).floor)).abs < 1e-4
	    #  n/4
		s += n.to_s + "/4"
	  elsif ((a * 6) - (n = (a * 6 + 0.5).floor)).abs < 1e-4
	    #  n/6
		s += n.to_s + "/6"
	  else
	    s += sprintf("%.4f", a)
	  end
	end
	["x", "y", "z"].each_with_index { |t, i|
	  a = sym[i, j]
	  if a.abs < 1e-4
	    next
	  elsif (a - 1.0).abs < 1e-4
	    s += (s == "" ? t : "+" + t)
	  elsif (a + 1.0).abs < 1e-4
	    s += "-" + t
	  else
	    if a > 0.0 && s != ""
		  s += "+"
		elsif a < 0.0
		  s += "-"
		end
		s += sprintf("%.4f", a.abs) + t
	  end
	}
	str += s
	if j < 2
	  str += ", "
	end
  }
  str
end

def string_to_symmetry(str)
  begin
    sary = str.split(/, */)
    raise if sary.count != 3
	a = []
	sary.each_with_index { |s, j|
	  if s[0] != "-"
	    s = "+" + s
	  end
	  while s.length > 0
	    raise if s !~ /^([-+][.0-9\/]*)([xyzXYZ]?)/
		sa = Regexp.last_match[1]
		st = Regexp.last_match[2]
		s = Regexp.last_match.post_match
	    case st
		when "x", "X"
		  i = 0
		when "y", "Y"
		  i = 1
		when "z", "Z"
		  i = 2
		else
		  i = 3
		end
		raise if a[i * 3 + j] != nil
		if sa == "-"
		  aa = -1
		elsif sa == "+"
		  aa = 1
		elsif sa.index("/")
		  sa0, sa1 = sa.split("/")
		  aa = sa0.to_f / sa1.to_f
		else
		  aa = sa.to_f
		end
		a[i * 3 + j] = aa
	  end
	}
	12.times { |i|
	  a[i] ||= 0.0
	}
	return Transform.new(a)
  rescue
    raise "Cannot convert to symmetry operation: #{str}"
  end
end

end

class Molecule

#  Export ortep to File fp
#  If attr is a hash, it represents options for drawing.
#  "atoms"=>[[group, type, color, rad], [group, type, color, rad], ...]
#    group is an IntGroup, representing a group of atoms.
#    type is an atom type: 0 = boundary, 1 = boundary + principal, 2 = 1 + axes, 3 = 2 + shades, 4 = fixed size
#    color is 0..7 (0,1=black, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 7=yellow)
#    rad is the radius of the atom (in Angstrom) whose type is fixed size
#    If an atom appears in multiple entries, the later entry is valid.
#  "bonds"=>[[group1, group2, type, color], [group1, group2, type, color], ...]
#    group1 and group2 are IntGroups, reprensenting the atoms constituting the bonds.
#    type is a bond type: 0 to 4 for bonds with no shades to 4 shades.
#    color is 0..7 as in atoms.
#    If a bond appears in multiple entries, the later entry is valid.
def export_ortep(fp, attr = nil)

  #  Create atom list
  hidden = atom_group { |ap| !is_atom_visible(ap.index) }
  hydrogen = self.show_hydrogens
  expanded = self.show_expanded
  atomlist = atom_group { |ap|
    (ap.element != "H" || hydrogen) &&
    (ap.symop == nil || expanded) &&
    (!hidden.include?(ap.index))
  }

  #  Title
  fp.printf "%-78.78s\n", self.name + ": generated by Molby at " + Time.now.to_s

  #  Cell parameters
  cp = self.cell
  if cp == nil
    cp = [1, 1, 1, 90, 90, 90]
  end
  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]
  
  #  Symmetry operations
  syms = self.symmetries
  if syms == nil || syms.length == 0
   fp.print "1             0  1  0  0              0  0  1  0              0  0  0  1\n"
  else
    syms.each_with_index { |s, i|
      a = s.to_a
      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]
    }
  end

  #  Atoms (all symmetry unique atoms regardless they are visible or not)
  n = 0
  aattr = (attr ? attr["atoms"] : nil)
  each_atom { |ap|
    break if ap.symop != nil
    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
	rad = -1.0
	if aattr
	  aattr.reverse_each { |at|
	    if at[0].include?(ap.index)
		  if at[1] == 4
		    rad = at[3].to_f
		    if rad == 0.0
		      rad = 0.1
		    end
		  end
		  break
		end
	  }
	end
    an = ap.aniso
    if an != nil && rad < 0.0
      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
    else
	  type = 7
	  if rad < 0.0
        rad = ap.temp_factor
        rad = 1.2 if rad <= 0
	    type = 6
	  end
      fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", rad, 0.0, 0.0, 0.0, 0.0, 0.0, type
    end
    n += 1
  }
  natoms_tep = n

  #  Special points to specify cartesian axes
  axis, angle = self.get_view_rotation
  tr = Transform.rotation(axis, -angle)
  org = self.get_view_center
  x = org + tr.column(0)
  y = org + tr.column(1)
  tr = self.cell_transform
  if tr
    tr = tr.inverse
  else
    tr = Transform.identity
  end
  org = tr * org
  x = tr * x
  y = tr * y
  fp.printf " CNTR                      %9.4f%9.4f%9.4f        0\n", org.x, org.y, org.z
  fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
  fp.printf " X                         %9.4f%9.4f%9.4f        0\n", x.x, x.y, x.z
  fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
  fp.printf " Y                         %9.4f%9.4f%9.4f        0\n", y.x, y.y, y.z
  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

  #  Initialize
  fp.print  "      201\n"

  #  Pen size
  fp.print  "      205        8\n"

  #  Paper size, margin, viewing distance
  fp.print  "      301      6.6      6.6        0      0.0\n"

  #  Sort atoms by symop and index
  acodes = Hash.new
  ig = IntGroup.new   #  Used for collecting contiguous atoms
  each_atom(atomlist) { |ap|
    idx, symcode = ap.to_adc
    acode = symcode * self.natoms + idx - 1
    acodes[acode] = ap.index
    ig.add(acode)
  }
  index2an = []       #  Index in Molecule to Index in ORTEP
  ig.each_with_index { |acode, i|
    index2an[acodes[acode]] = i + 1
  }
  i = 0
  adcs = []
  while (r = ig.range_at(i)) != nil
    s = r.first
    s = (s / self.natoms) + (s % self.natoms + 1) * 100000  #  Rebuild true ADC (atom index is in the upper digit)
    e = r.last
    e = (e / self.natoms) + (e % self.natoms + 1) * 100000
    if s < e
      adcs.push(s)
      adcs.push(-e)
    else
      adcs.push(s)
    end
    i += 1
  end
  k = 0

  #  Atom list
  adcs.each_with_index { |a, i|
    if k == 0
      fp.print "      401"
    end
    fp.printf "%9d", a
    k += 1
    if i == adcs.length - 1 || k == 6 || (k == 5 && i < adcs.length - 2 && adcs[i + 2] < 0)
      fp.print "\n"
      k = 0
    end
  }

  #  Axes
  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
#  fp.print  "      502        1      0.0        2      0.0        3      0.0\n"
  
  #  Autoscale
  fp.print  "      604                               1.538\n"
  
  #  Explicit bonds
  bond_inst = Array.new(35) { [] }   #  Bonds for types 0 to 4 and colors 1 to 7
  battr = (attr ? attr["bonds"] : nil)
  bonds.each { |b|
    next if !atomlist.include?(b[0]) || !atomlist.include?(b[1])
	btype = bcol = nil
	if battr
	  battr.reverse_each { |at|
	    if (at[0].include?(b[0]) && at[1].include?(b[1])) || (at[0].include?(b[1]) && at[1].include?(b[0]))
		  btype = at[2] % 5
		  bcol = (at[3] != 0 ? at[3] - 1 : at[3]) % 7
		  break
		end
	  }
	end
	if btype == nil
	  bcol = 0
      an1 = atoms[b[0]].atomic_number
      an2 = atoms[b[1]].atomic_number
      if an1 == 1 || an2 == 1
        btype = 0
      elsif an1 <= 8 && an2 <= 8
        btype = 2
      else
        btype = 4
      end
	end
    bond_inst[btype * 7 + bcol].push(b[0], b[1])
  }

  #  Output bond specifications
  #  Avoid including too many ADCs in a single 811/821 instruction
  #  (Upper limit is 140. Here we divide at every 36 ADCs)
  output_bonds = lambda { |icode|
    bond_inst.each_with_index { |inst, ii|
      next if inst.length == 0
	  btype = ii / 7 + 1
	  if icode / 10 == 81   #  811 instructions
	    fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
	  end
      inst.each_with_index { |b, i|
        if i % 6 == 0
          fp.printf "  %d   %3s", (i >= inst.length - 6 || i % 36 == 30 ? 2 : 1), (i % 36 == 0 ? icode.to_s : "")
        end
        idx, scode = atoms[b].to_adc
        fp.printf "%9d", idx * 100000 + scode
        if i % 6 == 5 || i == inst.length - 1
          fp.print "\n"
          if i == inst.length - 1 || i % 36 == 35
            fp.printf "%21s%3d%12s%6.3f\n", "", btype, "", 0.05
          end
        end
      }
    }
  }

  fp.print "  0  1001     0.02\n"  #  Activate hidden line removal
  output_bonds.call(821)
  
  #  Atom types
  #  Atom type 0=714 (boundary), 1=712 (+principal), 2=716 (+axes), 3=711 (+shades)
  atom_inst = Array.new(28) { IntGroup.new }
  aattr = (attr ? attr["atoms"] : nil)
  atomlist.each { |i|
    atype = acol = nil
	if aattr
	  aattr.reverse_each { |at|
	    if at[0].include?(i)
		  atype = at[1] % 4
		  acol = (at[2] != 0 ? at[2] - 1 : at[2]) % 7
		  break
		end
	  }
	end
	if atype == nil
	  acol = 0
      an1 = atoms[i].atomic_number
      if an1 == 1
        atype = 0
      elsif an1 <= 6
        atype = 1
      else
        atype = 3
      end
	end
    idx, scode = atoms[i].to_adc
    atom_inst[atype * 7 + acol].add(idx)
  }
  (atom_inst.count - 1).downto(0) { |ii|
    inst = atom_inst[ii]
	next if inst.count == 0
	fp.printf "      204%9d\n", ii % 7 + 1   #  Pen color
    i = 0
	atype = [714, 712, 716, 711][ii / 7]
    while (r = inst.range_at(i)) != nil
      fp.printf "  1   %3d\n", atype
      fp.printf "%27s%9d%9d\n", "", r.first, r.last
      i += 1
    end
  }

  output_bonds.call(811)

  #  Close plot
  fp.print "      202\n"
  fp.print "       -1\n"
  
end

def savetep(filename)
  if natoms == 0
	raise MolbyError, "cannot save ORTEP input; the molecule is empty"
  end
  fp = open(filename, "wb")
  export_ortep(fp)
  fp.close
  return true
end

end

#  Best-fit planes
#  Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
#       T. Ito, Acta Cryst 1981, A37, 621-624

#  An object to describe the best-fit plane
#  
#  Contains the plane coefficients and the constant (a, b, c, d for ax + by + cz + d = 0),
#  the error matrix, and the metric tensor.
#
class Molby::Plane
  attr_accessor :molecule, :group, :coeff, :const, :error_matrix, :metric_tensor
  def initialize(mol, group, coeff, const, err, met)
    @molecule = mol
    @group = group
    @coeff = coeff
    @const = const
    @error_matrix = err
    @metric_tensor = met
    self
  end
  def sigma
    [sqrt_safe(@error_matrix[0, 0]), sqrt_safe(@error_matrix[1, 1]), sqrt_safe(@error_matrix[2, 2]), sqrt_safe(@error_matrix[3, 3])]
  end
  def inspect
    s = sprintf("Molby::Plane[\n coeff, const = [[%f, %f, %f], %f],\n", @coeff.x, @coeff.y, @coeff.z, @const)
    s += sprintf(" sigma = [[%10.4e, %10.4e, %10.4e], %10.4e],\n", *self.sigma)
    (0..3).each { |i|
      s += (i == 0 ? " error_matrix = [" : "     ")
      (0..i).each { |j|
        s += sprintf("%12.6e%s", @error_matrix[j, i], (j == i ? (i == 3 ? "],\n" : ",\n") : ","))
      }
    }
    s += sprintf(" molecule = %s\n", @molecule.inspect)
    s += sprintf(" group = %s\n", @group.inspect)
    (0..3).each { |i|
      s += (i == 0 ? " metric_tensor = [" : "     ")
      (0..3).each { |j|
        s += sprintf("%12.6e%s", @metric_tensor[j, i], (j == 3 ? (i == 3 ? "]]\n" : ",\n") : ","))
      }
    }
    s
  end
  def distance(ap)
    if ap.is_a?(AtomRef)
      fr = ap.fract_r
      sig = ap.sigma
    else
      fr = Vector3D[*ap]
      sig = Vector3D[0, 0, 0]
    end
    d = fr.dot(@coeff) + @const
    sig1 = (@coeff.x * sig.x) ** 2 + (@coeff.y * sig.y) ** 2 + (@coeff.z * sig.z) ** 2
    sig2 = LAMatrix.multiply("t", fr, @error_matrix, fr)[0, 0]
    if ap.is_a?(AtomRef) && ap.molecule == @molecule && @group.include?(ap.index)
      #  The atom defines the plane
      sig0 = sig1 - sig2
      sig0 = 0.0 if sig0 < 0.0
    else
      sig0 = sig1 + sig2
    end  
    return d, sqrt_safe(sig0)
  end
  def dihedral(plane)
    e1 = @error_matrix.submatrix(0, 0, 3, 3)
    e2 = plane.error_matrix.submatrix(0, 0, 3, 3)
    m = @metric_tensor.submatrix(0, 0, 3, 3)
    c = plane.coeff
    cos_t = plane.coeff.dot(m * @coeff)
    if cos_t > 1.0
      cos_t = 1.0
    elsif cos_t < -1.0
      cos_t = -1.0
    end
    t = acos(cos_t)
    sig_t = (m * e1).trace + (m * e2).trace
    if sig_t < t * t
      w = 1.0 / sin(t)
      sig_t = w * w * (c.dot(LAMatrix.multiply(m, e1, m, c)) + @coeff.dot(LAMatrix.multiply(m, e2, m, @coeff)))
    end
    t *= 180.0 / PI
    sig_t = sqrt_safe(sig_t) * 180.0 / PI
    return t, sig_t
  end
end

class Molecule

#  Calculate best-fit plane for the given atoms
#  Return value: a Molby::Plane object

def plane(group)

  #  Number of atoms
  dim = group.length

  #  Positional parameters and standard deviations
  x = []
  sig = []
  sig_min = 1e10
  each_atom(group) { |ap|
    x.push(ap.fract_r)
    sig.push(ap.sigma)
    if (s = ap.sigma_x) > 0.0 && s < sig_min
      sig_min = s
    end
    if (s = ap.sigma_y) > 0.0 && s < sig_min
      sig_min = s
    end
    if (s = ap.sigma_z) > 0.0 && s < sig_min
      sig_min = s
    end
  }
  if sig_min == 1e10
    sig_min = 1e-12
  end
  sig.each { |s|
    s.x = sig_min if s.x < sig_min
    s.y = sig_min if s.y < sig_min
    s.z = sig_min if s.z < sig_min
  }

  #  The metric tensor of the reciprocal lattice
  #  g[j, i] = (ai*).dot(aj*), where ai* and aj* are the reciprocal axis vectors
  t = self.cell_transform
  if t.nil?
    t = Transform.identity
  end
  g2inv = LAMatrix[t]
  g2 = g2inv.inverse
  g2[3, 3] = 0.0
  g2inv[3, 3] = 0.0
  g = LAMatrix.multiply("t", g2, g2)

  #  The variance-covariance matrices of the atomic parameters
  #  mm[k][n] is a 3x3 matrix describing the correlation between the atoms k and n,
  #  and its components are defined as: sigma_k[i] * sigma_n[j] * corr[k, i, n, j],
  #  where corr(k, i, n, j) is the correlation coefficients between the atomic parameters
  #  k[i] and n[j].
  mm = Array.new(dim) { Array.new(dim) }
  zero = LAMatrix.zero(3, 3)
  dim.times { |k|
    dim.times { |n|
      mkn = LAMatrix.new(3, 3)
      if k == n
        3.times { |i|
          3.times { |j|
            if i == j
              mkn[j, i] = sig[k][i] * sig[n][j]
            else
              #  Inter-coordinate correlation should be implemented here
            end
          }
        }
      else
        #  Inter-atomic correlation should be implemented here
      end
      mm[k][n] = (mkn == zero ? zero : mkn)
    }
  }

  #  The variance-covariance matrix of the atom-plance distances
  #  m[j, i] = v.transpose * mm[i][j] * v, where v is the plane coefficient vector
  #  The inverse of m is the weight matrix
  m = LAMatrix.new(dim, dim)
  
  #  The matrix representation of the atomic coordinates
  #  y[j, i] = x[i][j] (for j = 0..2), -1 (for j = 3)
  #  y * LAMatrix[a, b, c, d] gives the atom-plane distances for each atom
  y = LAMatrix.new(4, dim)
  dim.times { |i|
    y[0, i] = x[i].x
    y[1, i] = x[i].y
    y[2, i] = x[i].z
    y[3, i] = 1.0
  }

  #  The coefficients to be determined
  n0 = LAMatrix[1, 1, 1, 0]
  v = LAMatrix[1, 1, 1]     #  The coefficient part

  iter = 0
  while iter < 20

    iter += 1

    #  Set zero to the "constant" part, and normalize the "coefficient" part
    n0[0, 3] = 0.0
    n0 = g2 * n0
    n0.multiply!(1.0 / n0.fnorm)
    n0 = g2inv * n0
    3.times { |i| v[0, i] = n0[0, i] }

    #  Build the variance-covariance matrix    
    dim.times { |i|
      dim.times { |j|
        m[j, i] = LAMatrix.multiply("t", v, mm[i][j], v)[0, 0]
      }
    }
    c = LAMatrix.multiply("t", y, "i", m, y)

    #  Invert c: only the inverse is used in the following, so c is inversed destructively
    cinv = c.inverse!
 
    if iter == 1

      #  Determine the tentative solution, which is given by the eigenvector of cinv * g
      #  for the largest eigenvalue
      evals, evecs = (cinv * g).eigenvalues
      4.times { |i| n0[0, i] = evecs[3, i] }

    else

      #  Convert the coefficient vector to the reciprocal space
      h = g * n0
      
      #  Determine multiplier
      #  In this implementation, the sign of delta-n is opposite from that used in
      #  the reference
      lam = 1.0 / (LAMatrix.multiply("t", h, cinv, h)[0, 0])
      
      #  Solve the linearized equation
      #  (Is the equation 21 in the reference really correct? Shouldn't it read
      #   B = 1 - lambda * C.inverse * H* ? )
      b = LAMatrix.multiply(lam, cinv, g)
      b.sub!(LAMatrix.identity(4))

      dn = b * n0
      n0 += dn

      break if dn[0, 0] ** 2 + dn[0, 1] ** 2 + dn[0, 2] ** 2 < 1e-9

    end
  end

  #  Error matrix = b * cinv * b.transpose
  em = LAMatrix.multiply(b, cinv, "t", b)
  coeff = Vector3D[n0[0, 0], n0[0, 1], n0[0, 2]]
  const = n0[0, 3]

  return Molby::Plane.new(self, group, coeff, const, em, g)

end

def cmd_plane
  plane_settings = @plane_settings || Hash.new
  mol = self
  h = Dialog.new("Best-Fit Planes: " + mol.name, "Close", nil) {
    refresh_proc = lambda { |it|
      n = it[:tag][/\d/].to_i
      g = plane_settings["group#{n}"]
      if g
        str = g.inspect.sub!("IntGroup[", "").sub!("]", "")
        set_value("group#{n}", str)
        if n == 1 || n == 2
          p = mol.plane(g) rescue p = nil
          plane_settings["plane#{n}"] = p
          if p
            coeff = p.coeff
            const = p.const
            sig = p.sigma
            aps = (n == 1 ? "" : "'")
            str = sprintf("a%s = %f(%f)\nb%s = %f(%f)\nc%s = %f(%f)\nd%s = %f(%f)",
                          aps, coeff.x, sig[0],
                          aps, coeff.y, sig[1],
                          aps, coeff.z, sig[2],
                          aps, const, sig[3])
            set_value("result#{n}", str)
          else
            set_value("result#{n}", "")
          end
          p1 = plane_settings["plane1"]
          p2 = plane_settings["plane2"]
          if p1 && p2
            t, sig = p1.dihedral(p2)
            str = sprintf("%f(%f)", t, sig)
            set_value("dihedral", str)
          else
            set_value("dihedral", "")
          end
        else
          p = plane_settings["plane1"]
          if p
            str = ""
            mol.each_atom(g) { |ap|
              d, sig = p.distance(ap)
              str += sprintf("%s  %f(%f)\n", ap.name, d, sig)
            }
            str.chomp!
          else
            str = ""
          end
          set_value("result#{n}", str)
        end
      else
        set_value("group#{n}", "")
        set_value("result#{n}", "")
      end
    }
    set_proc = lambda { |it|
      n = it[:tag][/\d/].to_i
      sel = mol.selection
      if sel.count > 0
        str = sel.inspect.sub!("IntGroup[", "").sub!("]", "")
        set_value("group#{n}", str)
        plane_settings["group#{n}"] = sel
      else
        plane_settings["group#{n}"] = nil
      end
      refresh_proc.call(it)
    }
    text_proc = lambda { |it|
      n = it[:tag][/\d/].to_i
      str = it[:value].gsub(/[^-.,0-9]/, "")  #  Remove unsane characters
      g = eval("IntGroup[#{str}]") rescue g = nil
      plane_settings["group#{n}"] = g
      refresh_proc.call(it)
    }
    layout(3,
      item(:text, :title=>"Plane 1 (ax + by + cz + d = 0)"),
      -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group1", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button1", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result1"),     
      item(:button, :title=>"Recalculate", :tag=>"refresh1", :action=>refresh_proc),
      item(:line),
      -1, -1,
      item(:text, :title=>"Plane 2 (a'x + b'y + c'z + d' = 0)"),
      -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group2", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button2", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result2"),
      item(:button, :title=>"Recalculate", :tag=>"refresh2", :action=>refresh_proc),
      item(:text, :title=>"Dihedral angle with Plane 1"), -1, -1,
      -1,
      item(:textfield, :width=>240, :height=>20, :tag=>"dihedral"), -1,
      item(:line),
      -1, -1,
      item(:text, :title=>"Distance from Plane 1"), -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group3", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button3", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result3"),
      item(:button, :title=>"Recalculate", :tag=>"refresh3", :action=>refresh_proc)
    )
    refresh_proc.call(item_with_tag("refresh1"))
    refresh_proc.call(item_with_tag("refresh2"))
    refresh_proc.call(item_with_tag("refresh3"))
	show
  }
  @plane_settings = plane_settings
end

#  Calculate bond length and angles with standard deviations
#  args can be a single IntGroup (calculate bonds and angles including those atoms)
#  or arrays of 2 or 3 integers (explicitly specifying bonds and angles)
def bond_angle_with_sigma(*args)

  if args.length >= 2 || (args.length == 1 && !args[0].is_a?(IntGroup))
    #  Bonds and angles are explicitly specified
    bonds = []
    angles = []
    args.each { |arg|
      if arg.length == 2 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer)
        bonds.push(arg)
      elsif arg.length == 3 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer) && arg[2].is_a?(Integer)
        angles.push(arg)
      else
        raise MolbyError, "Invalid argument #{arg.inspect}"
      end
    }
  else
    if args.length == 0
	  g = nil
	else
	  g = args[0]
	end
	bonds = self.bonds.select { |b|
	  (g == nil || (g.include?(b[0]) && g.include?(b[1]))) &&
	    (atoms[b[0]].symop == nil || atoms[b[1]].symop == nil)
	}
	angles = self.angles.select { |ang|
	  (g == nil || (g.include?(ang[0]) && g.include?(ang[1]) && g.include?(ang[2]))) &&
	    (atoms[ang[0]].symop == nil || atoms[ang[1]].symop == nil || atoms[ang[2]].symop == nil)
	}
  end

  #  A list of interatomic distance, its differential coefficients, and other
  #  useful quantities.
  #  The hash tag is a list [i, j] (i and j are the atom indices) or [i, j, k]
  #  (i, j, k are the atom indices, and r(ijk) is the distance between the atom i
  #  and the center point between the atoms j and k).
  #  The value is a list of following quantities:
  #    index 0: rij or r(ijk)
  #    index 1-9: d(rij)/d(xij), d(rij)/d(yij), d(rij)/d(zij),
  #      d(rij)/da, d(rij)/db, d(rij)/dc, d(rij)/d(alpha), d(rij)/d(beta), d(rij)/d(gamma)
  #    index 10: the list of the "base atom"
  #    index 11: the list of the transform matrices
  dlist = Hash.new

  #  A list of fractional coordinates and sigmas
  fract = []
  sigma = []
  each_atom { |ap|
    fract.push(ap.fract_r)
    sigma.push(ap.sigma)
  }

  #  A list of base atoms (for symmetry-related atoms) and transform matrices
  bases = []
  trans = []
  symcode = []
  trans_i = Transform.identity
  each_atom { |ap|
    sym = ap.symop
    bases.push(sym ? sym[4] : ap.index)
	if sym
	  tr = transform_for_symop(sym).transpose
      tr[3, 0] = tr[3, 1] = tr[3, 2] = 0.0
	else
	  tr = trans_i
	end
    trans.push(tr)
	symcode.push(sym ? sprintf("%d_%d%d%d", sym[0] + 1, sym[1] + 5, sym[2] + 5, sym[3] + 5) : ".")
  }

  #  Unit cell parameter
  cell = self.cell
  if cell == nil
    cell = [1, 1, 1, 90, 90, 90, 0, 0, 0, 0, 0, 0]
  elsif cell.length == 6
    cell.push(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)  #  No sigma information
  end
  # $a, $b, $c, $alpha, $beta, $gamma, $sig_a, $sig_b, $sig_c, $sig_alpha, $sig_beta, $sig_gamma = self.cell
  cos_a = cos(cell[3] * PI / 180.0)
  cos_b = cos(cell[4] * PI / 180.0)
  cos_c = cos(cell[5] * PI / 180.0)
  abc = cell[0] * cell[1] * cos_c
  bca = cell[1] * cell[2] * cos_a
  cab = cell[2] * cell[0] * cos_b
  aa = cell[0] * cell[0]
  bb = cell[1] * cell[1]
  cc = cell[2] * cell[2]

  get_dlist = lambda { |_i, _j, _k|
    if _k != nil
      if _j > _k
        _j, _k = _k, _j
      end
      _p = dlist[[_i, _j, _k]]
      return _p if _p != nil
      _p = (dlist[[_i, _j, _k]] = [])
      _vij = fract[_i] - (fract[_j] + fract[_k]) * 0.5
      _p[10] = [bases[_i], bases[_j], bases[_k]]
      _p[11] = [trans[_i], trans[_j], trans[_k]]
    else
      if _i > _j
        _i, _j = _j, _i
      end
      _p = dlist[[_i, _j]]
	  return _p if _p != nil
      _p = (dlist[[_i, _j]] = [])
      _vij = fract[_i] - fract[_j]
      _p[10] = [bases[_i], bases[_j]]
      _p[11] = [trans[_i], trans[_j]]
    end
    _xij = _vij.x
    _yij = _vij.y
    _zij = _vij.z
    _dij = sqrt_safe(aa * _xij * _xij + bb * _yij * _yij + cc * _zij * _zij + 2 * bca * _yij * _zij + 2 * cab * _zij * _xij + 2 * abc * _xij * _yij)
    _p[0] = _dij
    _p[1] = (aa * _xij + abc * _yij + cab * _zij) / _dij
    _p[2] = (bb * _yij + bca * _zij + abc * _xij) / _dij
    _p[3] = (cc * _zij + cab * _xij + bca * _yij) / _dij
    _p[4] = (cell[0] * _xij * _xij + cell[2] * cos_b * _zij * _xij + cell[1] * cos_c * _xij * _yij) / _dij
    _p[5] = (cell[1] * _yij * _yij + cell[0] * cos_c * _xij * _yij + cell[2] * cos_a * _yij * _zij) / _dij
    _p[6] = (cell[2] * _zij * _zij + cell[1] * cos_a * _yij * _zij + cell[0] * cos_b * _zij * _xij) / _dij
    _p[7] = (-cell[1] * cell[2] * sin(cell[3] * PI / 180.0) * _yij * _zij) * (PI / 180.0) / _dij
    _p[8] = (-cell[2] * cell[0] * sin(cell[4] * PI / 180.0) * _zij * _xij) * (PI / 180.0) / _dij
    _p[9] = (-cell[0] * cell[1] * sin(cell[5] * PI / 180.0) * _xij * _yij) * (PI / 180.0) / _dij
    return _p
  }

  diff_by_rn = lambda { |_dl, _n, _ijk|
    #  dl = dlist(i, j)
    #  return value: Vector3D[ d(rij)/d(xn), d(rij)/d(yn), d(rij)/d(zn) ]
    #  If ijk is true, then dl is dlist(i, j, k)
    _dv = Vector3D[_dl[1], _dl[2], _dl[3]]
    _c = Vector3D[0, 0, 0]
    if _dl[10][0] == _n
      _c += _dl[11][0] * _dv
    end
    if _ijk
      if _dl[10][1] == _n
        _c -= _dl[11][1] * _dv * 0.5
      end
      if _dl[10][2] == _n
        _c -= _dl[11][2] * _dv * 0.5
      end
    else
      if _dl[10][1] == _n
        _c -= _dl[11][1] * _dv
      end
    end
	return _c
  }

  notate_with_sigma = lambda { |_val, _sig|
    if _sig == 0.0
      return sprintf "%.4f", _val
    end
    _lg = log(_sig.abs / 1.95) / log(10.0)
    _n = -(_lg.floor)
    _val2 = (_val * (10 ** _n) + 0.5).floor * (0.1 ** _n)
    return sprintf "%.#{_n}f(%d)", _val2, (_sig * (10 ** _n) + 0.5).floor
  }

  results = []

  c = Vector3D[0, 0, 0]
  bonds.each { |b|
    i = b[0]
    j = b[1]
    if i > j
      i, j = j, i
    end
    p = get_dlist.call(i, j, nil)
    sig = 0.0
    p[10].uniq.each { |k|
      s = sigma[k]
      next unless s
      c = diff_by_rn.call(p, k, false)
      #  Test
      if nil
        apk = atoms[k]
        r = apk.fract_r
        d0 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0.00001, 0, 0]
        amend_by_symmetry(p[10])
        d1 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0, 0.00001, 0]
        amend_by_symmetry(p[10])
        d2 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0, 0, 0.00001]
        amend_by_symmetry(p[10])
        d3 = calc_bond(i, j)
        apk.fract_r = r
        amend_by_symmetry(p[10])
        printf " dr/dv[%d] cal %10.5g %10.5g %10.5g\n", k, c[0], c[1], c[2]
        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
      end
      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]
    }
    6.times { |n|
      sig += (p[4 + n] * cell[6 + n]) ** 2
    }
    sig = sqrt_safe(sig)
    results.push([[i, j], notate_with_sigma.call(p[0], sig), symcode[i], symcode[j], nil])
  }

  angles.each { |ang|
    i = ang[0]
    j = ang[1]
    k = ang[2]
    p0 = get_dlist.call(i, j, nil)
    p1 = get_dlist.call(j, k, nil)
    p2 = get_dlist.call(i, k, nil)
    p3 = get_dlist.call(j, i, k)
    t123 = acos_safe((p0[0] ** 2 + p1[0] ** 2 - p2[0] ** 2) / (2 * p0[0] * p1[0]))
    t124 = acos_safe((p0[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p0[0] * p3[0]))
    t324 = acos_safe((p1[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p1[0] * p3[0]))
    dtdr12 = -(p0[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p0[0] ** 2) * sin(t124))
    dtdr23 = -(p1[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p1[0] ** 2) * sin(t324))
    dtdr13 = p2[0] / (sin(t124) * 4 * p0[0] * p3[0]) + p2[0] / (sin(t324) * 4 * p1[0] * p3[0])
    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))
    pp = (p0[10] + p1[10] + p2[10] + p3[10]).uniq
    sig = 0.0
    pp.each { |n|
      s = sigma[n]
      next unless s
      c = Vector3D[0, 0, 0]
      c += diff_by_rn.call(p0, n, false) * (dtdr12 * 180.0 / PI)
      c += diff_by_rn.call(p1, n, false) * (dtdr23 * 180.0 / PI)
      c += diff_by_rn.call(p2, n, false) * (dtdr13 * 180.0 / PI)
      c += diff_by_rn.call(p3, n, true) * (dtdr24 * 180.0 / PI)
      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]
    }
    dd = dtdr12 * p0[4] + dtdr23 * p1[4] + dtdr13 * p2[4] + dtdr24 * p3[4]
    sig += dd * dd * cell[6] * cell[6]
    dd = dtdr12 * p0[5] + dtdr23 * p1[5] + dtdr13 * p2[5] + dtdr24 * p3[5]
    sig += dd * dd * cell[7] * cell[7]
    dd = dtdr12 * p0[6] + dtdr23 * p1[6] + dtdr13 * p2[6] + dtdr24 * p3[6]
    sig += dd * dd * cell[8] * cell[8]
    dd = dtdr12 * p0[7] + dtdr23 * p1[7] + dtdr13 * p2[7] + dtdr24 * p3[7]
    sig += dd * dd * cell[9] * cell[9]
    dd = dtdr12 * p0[8] + dtdr23 * p1[8] + dtdr13 * p2[8] + dtdr24 * p3[8]
    sig += dd * dd * cell[10] * cell[10]
    dd = dtdr12 * p0[9] + dtdr23 * p1[9] + dtdr13 * p2[9] + dtdr24 * p3[9]
    sig += dd * dd * cell[11] * cell[11]
    sig = sqrt_safe(sig)
    results.push([[i, j, k], notate_with_sigma.call(t123*180/PI, sig), symcode[i], symcode[j], symcode[k]])
  }
  results
end

def cmd_bond_angle_with_sigma
  mol = self
  Dialog.new("Bond & Angle with Sigma:" + self.name, nil, nil) {
    values = []
	clicked = []
	sel = mol.selection
	on_add_bond = lambda { |it|
	  values.push([nil, nil, "-", nil, nil, nil, "-", 0])
	  item_with_tag("table")[:refresh] = true
	  item_with_tag("table")[:selection] = IntGroup[values.count - 1]
	  item_with_tag("dump")[:enabled] = true
	}
	on_add_angle = lambda { |it|
	  values.push([nil, nil, nil, nil, nil, nil, nil, 0])
	  item_with_tag("table")[:refresh] = true
	  item_with_tag("table")[:selection] = IntGroup[values.count - 1]
	  item_with_tag("dump")[:enabled] = true
	}
	on_get_value = lambda { |it, row, col|
	  v = values[row][col]
	  if col < 3
	    if v.is_a?(Integer)
	      v = mol.atoms[v].name
		elsif v == nil
		  v = "(select)"
		end
	  elsif (col >= 4 && col <= 6) && (vv = values[row][col - 4]).is_a?(Integer)
	    if (s = mol.atoms[vv].symop) != nil
		  v = sprintf("%d_%d%d%d", s[0] + 1, s[1] + 5, s[2] + 5, s[3] + 5)
		else
		  v = "."
		end
	  end
	  return v
	}
	on_selection_changed = lambda { |it|
	  s = it[:selection]
	  if s && s.count > 0
	    set_attr("delete", :enabled=>true)
	    set_attr("dump", :enabled=>true)
	  else
	    set_attr("delete", :enabled=>false)
		set_attr("dump", :enabled=>false)
	  end
	}
	on_delete = lambda { |it|
	  s = attr("table", :selection)
	  if s
	    s.reverse_each { |idx|
		  values.delete_at(idx)
		}
	  end
	  tbl = item_with_tag("table")
	  tbl[:selection] = IntGroup[]
	  tbl[:refresh] = true
	  on_selection_changed.call(tbl)
	}
	atom_name = lambda { |ap|
	  (ap.molecule.nresidues >= 2 ? "#{ap.res_seq}:" : "") + ap.name
	}
    layout(1,
	  item(:table, :width=>500, :height=>300, :tag=>"table",
	    :columns=>[["atom1", 60], ["atom2", 60], ["atom3", 60],
		  "value(sigma)", "symop1", "symop2", "symop3"],
		:on_count=> lambda { |it| values.count },
		:on_get_value=>on_get_value,
		:on_selection_changed=>on_selection_changed),
	  item(:text, :title=>"(1) Hit 'New Bond' or 'New Angle', and (2) click on the atoms to measure in the main window.",
	    :font=>[10]),
	  layout(3,
	    item(:button, :title=>"New Bond", :width=>80, :font=>[10], :action=>on_add_bond),
		item(:button, :title=>"New Angle", :width=>80, :font=>[10], :action=>on_add_angle),
		item(:button, :title=>"Delete", :width=>80, :font=>[10], :action=>on_delete,
		  :tag=>"delete", :enabled=>false),
		:padding=>5, :margin=>0),
	  layout(2,
	    item(:view, :width=>480, :height=>1), -1,
	    item(:button, :title=>"Export to Clipboard", :tag=>"dump",
		  :action=>lambda { |item|
		    s = ""
			sel = attr("table", :selection)
			return if sel == nil || sel.count == 0
			sel.each { |j|
			  ss = ""
			  7.times { |i|
			    sss = on_get_value.call(item_with_tag("table"), j, i).to_s
				ss += (sss == "-" ? "" : sss) + (i == 6 ? "\n" : "  ")
			  }
			  s += ss
			}
			export_to_clipboard(s)
		  },
		  :enabled=>false),
		[item(:button, :title=>"Close", :action=>lambda { |item| hide }), {:align=>:right}])
	)
	on_document_modified = lambda { |*d|
	  newsel = mol.selection - sel
	  sel = mol.selection
	  row = (item_with_tag("table")[:selection] || [nil])[-1]
	  return unless row
	  if newsel.count == 1
	    #  Add new atom to the selected row
		val = newsel[0]
		lim = (values[row][2] == "-" ? 2 : 3)
		i = 0
		while i < lim
		  if values[row][i] == nil
		    values[row][i] = val
			break
		  end
		  i += 1
		end
		if i == lim
		  (lim - 1).times { |j|
		    values[row][j] = values[row][j + 1]
		  }
		  values[row][lim - 1] = val
		end
		if i < lim - 1
		  val1 = nil
		else
		  val1 = mol.bond_angle_with_sigma(values[row][0...lim])[0][1]
		end
		values[row][3] = val1
		item_with_tag("table")[:refresh] = true
	  end
	}
    listen(mol, "documentModified", on_document_modified)
	listen(mol, "documentWillClose", lambda { |*d| hide } )
	on_document_modified.call
	show
  }
end

def find_expanded_atom(base, symop)
  return base unless symop
  symopb = symop + [base]
  self.each_atom { |ap|
    if ap.symop == symopb
      return ap.index
    end
  }
  nil
end

def complete_by_symmetry
  if self.box == nil
    raise "Unit cell should be given"
  end
  verbose = false
  avec, bvec, cvec = self.box
  syms = []
  self.nsymmetries.times { |n|
    syms.push(transform_for_symop([n, 0, 0, 0], true))
  }
  frags = []
  self.each_fragment { |f|
    frags.push(f)
  }
  close_pairs = []
  self.each_atom { |ap|
    #  Find if ap is contained in expansion
	next if ap.symop != nil  #  Already expanded
    rad = Parameter.builtin.elements[ap.atomic_number].radius
    syms.each_with_index { |sym, sym_idx|
      27.times { |n|
        dx = n % 3 - 1
        dy = (n / 3) % 3 - 1
        dz = (n / 9) % 3 - 1
        next if dx == 0 && dy == 0 && dz == 0 && sym_idx == 0
        symop = [sym_idx, dx, dy, dz]
		next if self.find_expanded_atom(ap.index, symop)
        r = sym * ap.r + avec * dx + bvec * dy + cvec * dz
        close_atoms = self.find_close_atoms(r, 1.2, rad)
        if close_atoms.length > 0
          #  ap * [sym, dx, dy, dz] is included in the expansion
          close_atoms.each { |ca|
            next if ca > ap.index
            i1 = frags.index { |f| f.include?(ca) }
            i2 = frags.index { |f| f.include?(ap.index) }
            pp = close_pairs.find { |p1| p1[0] == i1 && p1[1] == i2 && p1[2] == symop }
			if pp == nil
			  pp = [i1, i2, symop, [], [], [], []]
			  close_pairs.push(pp)
		    end
            if (r - atoms[ca].r).length2 > 0.0001
			  #  Normal case (bond between ca and ap)
			  pp[3].push(ca)
			  pp[4].push(ap.index)
			else
              #  Special position
              pp[5].push(ca)
              pp[6].push(ap.index)
            end
          }
        end
      }
    }
  }
  puts "close_pairs = #{close_pairs}" if verbose
  expand_pairs = lambda { |i1, symop, depth|
    #  Find expanded fragment that is close to fragment [i1, symop]
    next [] if depth > 16
    retval = []
    close_pairs.each { |pp|
      next if i1 != nil && pp[0] != i1
	  next if pp[3].count == 0   #  No real expansion
      #  Multiply two symops
      if symop
        symop2 = symop_for_transform(transform_for_symop(pp[2]) * transform_for_symop(symop))
        puts "multiply two symops: #{pp[2].inspect} * #{symop.inspect} -> #{symop2.inspect}" if verbose
      else
        symop2 = pp[2]
      end
      if symop2[0] == 0
        #  Translation only
        if symop2[1] == 0 && symop2[2] == 0 && symop2[3] == 0
		  #  No expansion, but pp[3] and pp[4] are still needed to make bonds
		  symop2 = nil
	    end
		#  Expand the atoms only at the marginal
	    retval.push([nil, symop, symop2, pp[3], pp[4], pp[5], pp[6]])
      else
        #  Expand fragment
        retval.push([pp[1], symop, symop2, pp[3], pp[4], pp[5], pp[6]])
        retval.concat(expand_pairs.call(pp[1], symop2, depth + 1))
      end
    }
    next retval
  }
  ex_pairs = expand_pairs.call(nil, nil, 0)
  puts "ex_pairs = #{ex_pairs}" if verbose

  #  Expand fragments
  duplicates = []
  ex_pairs.each { |ex|
    #  ex[0]: fragment to expand, ex[1], ex[2]: symop,
    #  ex[3], ex[4]: atoms to make bonds, ex[5], ex[6]: duplicate atoms
    n1 = self.natoms
    if ex[0] == nil
      f = IntGroup[ex[4]]
    else
      f = frags[ex[0]]
    end
	if ex[2] != nil
      self.expand_by_symmetry(f, ex[2][0], ex[2][1], ex[2][2], ex[2][3], true)
      puts "expand(#{f}, #{ex[2]}) -> atoms[#{n1}..#{self.natoms - n1}]" if verbose
	end
    ex[3].each_with_index { |x, i|
      x1 = self.find_expanded_atom(x, ex[1])
	  if ex[2] != nil
        x2 = f.index(ex[4][i]) + n1  #  New index of the expanded atom
	  else
	    x2 = ex[4][i]   #  Base atom
	  end
      create_bond(x1, x2)
      puts "create_bond(#{x1}, #{x2})" if verbose
    }
    #  Register duplicate atoms
	if ex[2] != nil
      ex[5].each_with_index { |x, i|
        x1 = self.find_expanded_atom(x, ex[1])
        x2 = f.index(ex[6][i]) + n1
        duplicates.each { |d|
          if d.include?(x1)
            d.push(x2)
            x2 = nil
            break
          end
        }
        if x2
          duplicates.push([x1, x2])
        end
      }
    end
  }
  puts "duplicates = #{duplicates}" if verbose
  
  #  Remove duplicate atoms
  rem = []
  duplicates.each { |d|
    d.each_with_index { |dn, i|
      next if i == 0
      #  Remake bonds
      self.atoms[dn].connects.each { |nn|
        create_bond(d[0], nn)
        puts "create_bond(#{d[0]}, #{nn})" if verbose
      }
      rem.push(dn)
    }
  }
  remove(rem)
  puts "remove(#{rem})" if verbose
  
end

def create_packing_diagram
  if self.box == nil
    error_message_box "Unit cell is not defined."
	return
  end
  expansion_box = (@expansion_box ||= [0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
  h = Dialog.run("Create Packing Diagram") {
    layout(4,
	  item(:text, :title=>"Specify the expansion limits for each axis:"),
	  -1, -1, -1,
	  item(:text, :title=>"x"),
	  item(:textfield, :width=>80, :tag=>"xmin", :value=>sprintf("%.1f", expansion_box[0].to_f)),
	  item(:text, :title=>"..."),
	  item(:textfield, :width=>80, :tag=>"xmax", :value=>sprintf("%.1f", expansion_box[1].to_f)),
	  item(:text, :title=>"y"),
	  item(:textfield, :width=>80, :tag=>"ymin", :value=>sprintf("%.1f", expansion_box[2].to_f)),
	  item(:text, :title=>"..."),
	  item(:textfield, :width=>80, :tag=>"ymax", :value=>sprintf("%.1f", expansion_box[3].to_f)),
	  item(:text, :title=>"z"),
	  item(:textfield, :width=>80, :tag=>"zmin", :value=>sprintf("%.1f", expansion_box[4].to_f)),
	  item(:text, :title=>"..."),
	  item(:textfield, :width=>80, :tag=>"zmax", :value=>sprintf("%.1f", expansion_box[5].to_f)))
  }
  if h[:status] == 0
    @expansion_box[0] = h["xmin"].to_f
    @expansion_box[1] = h["xmax"].to_f
    @expansion_box[2] = h["ymin"].to_f
    @expansion_box[3] = h["ymax"].to_f
    @expansion_box[4] = h["zmin"].to_f
    @expansion_box[5] = h["zmax"].to_f
  else
    return
  end
  
  #  Enumerate all atoms within the expansion box
  syms = self.symmetries
  h = Hash.new
  xmin = @expansion_box[0]
  xmax = @expansion_box[1]
  ymin = @expansion_box[2]
  ymax = @expansion_box[3]
  zmin = @expansion_box[4]
  zmax = @expansion_box[5]
  xmin_d = xmin.floor
  xmax_d = xmax.floor - 1
  ymin_d = ymin.floor
  ymax_d = ymax.floor - 1
  zmin_d = zmin.floor
  zmax_d = zmax.floor - 1
  syms.each_with_index { |sym, i|
    each_atom { |ap|
      fr = sym * ap.fract_r
	  dx = fr.x.floor
	  dy = fr.y.floor
	  dz = fr.z.floor
	  fr.x -= dx
	  fr.y -= dy
	  fr.z -= dz
	  symopi = i * 1000000 + (50 - dx) * 10000 + (50 - dy) * 100 + 50 - dz
	  (h[symopi] ||= []).push(ap.index)
	  xmin_d.upto(xmax_d) { |xd|
	    next if fr.x + xd < xmin || fr.x + xd > xmax
	    ymin_d.upto(ymax_d) { |yd|
	      next if fr.y + yd < ymin || fr.y + yd > ymax
		  zmin_d.upto(zmax_d) { |zd|
		    next if fr.z + zd < zmin || fr.z + zd > zmax
			next if xd == 0 && yd == 0 && zd == 0
		    symopid = symopi + xd * 10000 + yd * 100 + zd
			(h[symopid] ||= []).push(ap.index)
		  }
		}
	  }
	}
  }
  h.keys.sort.each { |key|
    sym = key / 1000000
	dx = key % 1000000 / 10000 - 50
	dy = key % 10000 / 100 - 50
	dz = key % 100 - 50
	next if sym == 0 && dx == 0 && dy == 0 && dz == 0
#	puts "[#{sym},#{dx},#{dy},#{dz}]: #{h[key].inspect}"
	expand_by_symmetry(IntGroup[h[key]], sym, dx, dy, dz)
  }
end

def cmd_show_ortep
  mol = self
  begin
    tmp = create_temp_dir("ortep", mol.name)
  rescue
    error_message_box($!.to_s)
	return
  end
  tepexe = "#{ResourcePath}/ortep3/ortep3"
  if $platform == "win"
    tepexe += ".exe"
  end
  tepdata = []
  tepbounds = [0, 0, 400, 400]
  descs = {
	"Atoms"=>[
	  ["all", "Shaded", "black", ""],
	  ["Li-Ar", "Principals", "black", ""],
	  ["C", "Boundary", "black", ""],
	  ["H", "Fixed", "black", "0.1"]
	],
	"Bonds"=>[
	  ["all", "all", "4 Shades", "black"],
	  ["Li-Ar", "Li-Ar", "3 Shades", "black"],
	  ["C", "Li-Ar", "2 Shades", "black"],
	  ["H", "all", "Open", "black"]
	]
  }
  Dialog.new("Show ORTEP:" + mol.name, nil, nil, :resizable=>true, :has_close_box=>true) {
    tepview = nil  #  Forward declaration
	tab = "Atoms"
	columns = {
	  "Atoms"=>[["atom list", 65], ["type", 80], ["color", 40], ["radius", 50]],
	  "Bonds"=>[["list1", 65], ["list2", 65], ["type", 80], ["color", 40]]
	}
	atom_types = ["Boundary", "Principals", "Axes", "Shaded", "Fixed"]
	bond_types = ["Open", "1 Shade", "2 Shades", "3 Shades", "4 Shades"]
	colors = ["black", "red", "green", "blue", "cyan", "magenta", "yellow"]
	#  ORTEP attributes
	get_ortep_attr = lambda {
	  attr = Hash.new
	  make_atom_list = lambda { |s, ln|
	    ss = s.scan(/[-.0-9A-Za-z]+/)
		# puts ss.inspect
		ss.inject(IntGroup[]) { |r, it|
		  if it == "all"
		    r.add(mol.all)
			next r
		  elsif it =~ /-|(\.\.)/
		    s0 = Regexp.last_match.pre_match
			s1 = Regexp.last_match.post_match
		  else
		    s0 = s1 = it
		  end
		  if s0 =~ /^[0-9]+$/
		    #  Atomic numbers
			s0 = s0.to_i
			if s1 !~ /^[0-9]+$/ || (s1 = s1.to_i) < s0
			  raise "Bad atom list specification: #{s} in line #{ln}"
			end
			r.add(s0..s1)
		  else
		    #  Atomic symbols
			s0 = Parameter.builtin.elements.find_index { |p| p.name == s0.capitalize }
			s1 = Parameter.builtin.elements.find_index { |p| p.name == s1.capitalize }
			if s0 == nil || s1 == nil
			  raise "Bad element symbol specification: #{s} in line #{ln}"
			end
			(s0..s1).each { |an|
			  r.add(mol.atom_group { |ap| ap.atomic_number == an } )
			}
		  end
		  r
		}
      }
	  #  Atoms
	  at = []
	  descs["Atoms"].each_with_index { |d, idx|
		list = make_atom_list.call(d[0], idx)
		type = atom_types.find_index { |it| it == d[1] } || 0
		col = (colors.find_index { |it| it == d[2] } || 0) + 1
		rad = d[3]
		at.push([list, type, col, rad])
	  }
	  attr["atoms"] = at
	  #  Bonds
	  bd = []
	  descs["Bonds"].each_with_index { |d, idx|
		list1 = make_atom_list.call(d[0], idx)
		list2 = make_atom_list.call(d[1], idx)
		type = bond_types.find_index { |it| it == d[2] } || 0
		col = (colors.find_index { |it| it == d[3] } || 0) + 1
		bd.push([list1, list2, type, col])
	  }
	  attr["bonds"] = bd
	  attr
	}
    on_update_ortep = lambda { |it|
	  #  ORTEP attributes
	  attr = get_ortep_attr.call
	  #  Create ORTEP input in the temporary directory
	  open(tmp + "/TEP.IN", "w") { |fp| mol.export_ortep(fp, attr) }
	  #  Run ORTEP
	  cwd = Dir.pwd
	  Dir.chdir(tmp)
	  if FileTest.exist?("TEP001.PRN")
	    File.unlink("TEP001.PRN")
	  end
	  pid = call_subprocess(tepexe, "Running ORTEP", nil, "NUL", "NUL")
	  if FileTest.exist?("TEP001.PRN")
	    File.rename("TEP001.PRN", "TEP001.ps")
	  end
	  Dir.chdir(cwd)
	  if pid != 0
	    msg = "ORTEP execution in #{tmp} failed with status #{pid}."
	    message_box(msg, "ORTEP Failed", :ok, :warning)
	  else
	    open(tmp + "/TEP001.ps", "r") { |fp|
		  tepdata.clear
		  tepbounds = [100000, 100000, -100000, -100000]
		  fp.each_line { |ln|
		    ln.chomp!
		    x, y, c = ln.split
			if ln =~ /setrgbcolor/
			  tepdata.push(["color", x.to_f, y.to_f, c.to_f])
			  next
			elsif ln =~ /setgray/
			  x = x.to_f
			  tepdata.push(["color", x, x, x])
			  next
			end
			x = x.to_i
			y = y.to_i
			if c == "m"
			  tepdata.push([x, y])
			elsif c == "l"
			  tepdata[-1].push(x, y)
			else
			  next
			end
			tepbounds[0] = x if x < tepbounds[0]
			tepbounds[1] = y if y < tepbounds[1]
			tepbounds[2] = x if x > tepbounds[2]
			tepbounds[3] = y if y > tepbounds[3]
		  }
		  #  [x, y, width, height]
		  tepbounds[2] -= tepbounds[0]
		  tepbounds[3] -= tepbounds[1]
		}
		tepview.refresh_rect(tepview[:frame])
	  end
	}
	on_draw_ortep = lambda { |it|
	  clear
	  frame = it[:frame]
	  brush(:color=>[1, 1, 1, 1])
	  pen(:color=>[1, 1, 1, 1])
	  draw_rectangle(frame)
	  pen(:color=>[0, 0, 0, 1])
	  rx = (frame[2] - 10.0) / tepbounds[2]
	  ry = (frame[3] - 10.0) / tepbounds[3]
	  rx = ry if rx > ry
	  dx = (frame[2] - tepbounds[2] * rx) * 0.5
	  dy = (frame[3] + tepbounds[3] * rx) * 0.5
	  tepdata.each { |d|
	    if d[0] == "color"
		  pen(:color=>[d[1], d[2], d[3], 1])
		  next
		end
	    x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
		y0 = (dy - (d[1] - tepbounds[1]) * rx).to_i
		(d.length / 2 - 1).times { |i|
		  x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
		  y1 = (dy - (d[i * 2 + 3] - tepbounds[1]) * rx).to_i
		  draw_line(x0, y0, x1, y1)
		  x0 = x1
		  y0 = y1
		}
	  }
	}
	on_export_eps = lambda { |fname|
	  frame = item_with_tag("ortep")[:frame].dup
	  dx = dy = 5
	  rx = (frame[2] - dx * 2) / tepbounds[2]
	  ry = (frame[3] - dy * 2) / tepbounds[3]
	  rx = ry if rx > ry
	  frame[2] = (rx * tepbounds[2] + dx * 2).to_i
	  frame[3] = (rx * tepbounds[3] + dy * 2).to_i
	  #  Assumes A4 paper and center the bounding box
	  frame[0] = (595 - frame[2]) / 2
	  frame[1] = (842 - frame[3]) / 2
	  xmax = frame[0] + frame[2]
	  ymax = frame[1] + frame[3]
	  open(fname, "w") { |fp|
	    fp.print "%!PS-Adobe-3.0 EPSF-3.0
%%Creator: Molby
%%BoundingBox: #{frame[0]} #{frame[1]} #{xmax} #{ymax}
%%Pages: 1
%%Orientation: Landscape
%%BeginProlog
/m {moveto} def
/l {lineto} def
%%EndProlog
%%Page: 1 1
%%BeginPageSetup
0 setgray 1 setlinecap 1 setlinewidth
%%EndPageSetup
"
	    tepdata.each { |d|
	      if d[0] == "color"
		    fp.printf "%.3f %.3f %.3f setrgbcolor\n", d[1], d[2], d[3]
		    next
		  end
	      x0 = frame[0] + dx + (d[0] - tepbounds[0]) * rx
		  y0 = frame[1] + dy + (d[1] - tepbounds[1]) * rx
		  fp.printf "%.2f %.2f m\n", x0, y0
		  (d.length / 2 - 1).times { |i|
		    x1 = frame[0] + dx + (d[i * 2 + 2] - tepbounds[0]) * rx
		    y1 = frame[1] + dy + (d[i * 2 + 3] - tepbounds[1]) * rx
		    fp.printf "%.2f %.2f l\n", x1, y1
		    x0 = x1
		    y0 = y1
		  }
		  fp.print "stroke\n"
	    }
        fp.print "showpage\n%%Trailer\n%%EOF\n"
	  }
	}
	on_export_bitmap = lambda { |fname|
	  frame = item_with_tag("ortep")[:frame].dup
	  dx = dy = 5
	  scale = 5.0
	  rx = (frame[2] - dx * 2) * scale / tepbounds[2]
	  ry = (frame[3] - dy * 2) * scale / tepbounds[3]
	  rx = ry if rx > ry
	  frame[2] = (rx * tepbounds[2] + dx * 2).to_i
	  frame[3] = (rx * tepbounds[3] + dy * 2).to_i
      bmp = Bitmap.new(frame[2], frame[3], 32)
	  bmp.focus_exec {
	    clear
	    pen(:color=>[0, 0, 0, 1], :width=>scale)
	    tepdata.each { |d|
	      if d[0] == "color"
		    pen(:color=>[d[1], d[2], d[3], 1])
		    next
		  end
	      x0 = (dx + (d[0] - tepbounds[0]) * rx).to_i
		  y0 = (dy + (tepbounds[3] - d[1] + tepbounds[1]) * rx).to_i
		  (d.length / 2 - 1).times { |i|
		    x1 = (dx + (d[i * 2 + 2] - tepbounds[0]) * rx).to_i
		    y1 = (dy + (tepbounds[3] - d[i * 2 + 3] + tepbounds[1]) * rx).to_i
		    draw_line(x0, y0, x1, y1)
		    x0 = x1
		    y0 = y1
		  }
	    }
	  }
	  bmp.save_to_file(fname)
	}
	on_export = lambda { |it|
	  basename = (mol.path ? File.basename(mol.path, ".*") : mol.name)
      fname = Dialog.save_panel("Export ORTEP file:", mol.dir, basename + ".eps",
	    "Encapsulated PostScript (*.eps)|*.eps|PNG (*.png)|*.png|TIFF (*.tif)|*.tif|ORTEP Instruction (*.tep)|*.tep|All files|*.*")
	  return if !fname
	  ext = File.extname(fname).downcase
	  if ext == ".eps"
	    on_export_eps.call(fname)
	  elsif ext == ".png" || ext == ".tif" || ext == ".tiff"
	    on_export_bitmap.call(fname)
	  elsif ext == ".tep"
	    filecopy("#{tmp}/TEP.IN", fname)
      end
	}
	on_document_modified = lambda { |m|
	}
	#  Close handler (called when the close box is pressed or the document is closed)
	@on_close = lambda { |*d|
	  erase_old_logs(tmp)
	  tmp = nil
	  true
	}
	on_select_tab = lambda { |it|
	  tab = ["Atoms", "Bonds"][it[:value]]
	  # puts "tab = #{tab}"
	  table = item_with_tag("table")
	  table[:columns] = columns[tab]
	  table[:selection] = IntGroup[]
	  table[:refresh] = true
	}
	get_count = lambda { |it| descs[tab].count }
	get_value = lambda { |it, row, col| descs[tab][row][col] }
	has_popup_menu = lambda { |it, row, col|
	  if tab == "Atoms"
	    if col == 1
	      return atom_types
		elsif col == 2
		  return colors
		end
	  elsif tab == "Bonds"
		if col == 2
	      return bond_types
		elsif col == 3
		  return colors
		end
	  end
	  return nil
	}
	on_selection_changed = lambda { |it|
	  sel = it[:selection]
	  item_with_tag("remove")[:enabled] = (sel == nil || sel.count == 0 ? false : true)
	}
	is_item_editable = lambda { |it, row, col|
	  return has_popup_menu.call(it, row, col) == nil
	}
	popup_menu_selected = lambda { |it, row, col, sel|
	  titles = has_popup_menu.call(it, row, col)
	  return nil if titles == nil
	  descs[tab][row][col] = titles[sel]
	  it[:refresh] = true
	}
	set_value = lambda { |it, row, col, val|
	  descs[tab][row][col] = val
	  it[:refresh] = true	  
	}
	add_to_table = lambda { |it|
	  table = item_with_tag("table")
	  d = descs[tab][-1].dup
	  descs[tab].push(d)
	  table[:refresh] = true
	  table[:selection] = IntGroup[descs[tab].count - 1]
	  # puts descs[tab].inspect
	}
	remove_from_table = lambda { |it|
	  table = item_with_tag("table")
	  sel = table[:selection]
	  sel.reverse_each { |n|
	    descs[tab][n, 1] = []
	  }
	  table[:refresh] = true
	}
	layout(2,
	  layout(1,
	    item(:popup, :subitems=>["Atoms", "Bonds"], :action=>on_select_tab),
		item(:table, :width=>240, :height=>380, :flex=>[0,0,0,0,1,1], :tag=>"table",
		  :columns=>columns["Atoms"],
		  :on_count=>get_count,
		  :on_get_value=>get_value,
		  :on_set_value=>set_value,
		  :on_selection_changed=>on_selection_changed,
		  :is_item_editable=>is_item_editable,
		  :has_popup_menu=>has_popup_menu,
		  :on_popup_menu_selected=>popup_menu_selected),
		layout(2,
		  item(:button, :title=>"Add", :tag=>"add", :action=>add_to_table),
		  item(:button, :title=>"Remove", :tag=>"remove", :action=>remove_from_table, :enabled=>false),
		  :flex=>[0,1,0,0,0,0]),
	    :flex=>[0,0,0,0,0,1]),
	  layout(1,
	    item(:view, :width=>400, :height=>400, :tag=>"ortep", :on_paint=>on_draw_ortep, :flex=>[0,0,0,0,1,1]),
	    layout(2,
	      item(:button, :title=>"Refresh", :action=>on_update_ortep, :align=>:left),
		  item(:button, :title=>"Export as...", :action=>on_export, :align=>:right),
		  :flex=>[0,1,0,0,0,0]),
	    :flex=>[0,0,0,0,1,1]),
	  :flex=>[0,0,0,0,1,1]
	)
	set_min_size(480, 200)
	tepview = item_with_tag("ortep")
	listen(mol, "documentModified", on_document_modified)
	listen(mol, "documentWillClose", lambda { |m| close } )
	on_document_modified.call(nil)
    on_update_ortep.call(nil)
	show
  }
end

  def cmd_unit_cell
    mol = self
	new_box = nil
    hash = Dialog.run("Unit Cell") {
	  @mol = mol
	  @box = @mol.box
	  @box = (@box ? @box.dup : nil)
	  @box_save = nil
	  
	  def set_cell_value(item1)
		alpha = value("alpha").to_f * PI / 180.0
		beta = value("beta").to_f * PI / 180.0
		gamma = value("gamma").to_f * PI / 180.0
		if alpha.abs < 1e-8 || beta.abs < 1e-8 || gamma.abs < 1e-8
		  return
		end
	    if @box == nil
		  @box = [Vector3D[1,0,0],Vector3D[0,1,0],Vector3D[0,0,1],Vector3D[0,0,0],[1,1,1]]
		end
		av = Vector3D[1, 0, 0]
		bv = Vector3D[cos(gamma), sin(gamma), 0]
		cx = cos(beta)
		cy = (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma)
		cz = sqrt_safe(1.0 - cx * cx - cy * cy)
		cv = Vector3D[cx, cy, cz]
		x0 = @box[0].normalize
		z0 = (@box[0].cross(@box[1])).normalize
		y0 = (z0.cross(x0)).normalize
		tr = Transform[x0, y0, z0, Vector3D[0, 0, 0]]
		av = (tr * av) * value("a").to_f
		bv = (tr * bv) * value("b").to_f
		cv = (tr * cv) * value("c").to_f
		@box[0] = av
		@box[1] = bv
		@box[2] = cv
		3.times { |i|
		  3.times { |j|
		    if @box[i][j].abs < 1e-8
			  @box[i][j] = 0.0
			end
		  }
		}
		update_values(item1)
	  end
	  
	  def set_box_value(item1)
	    h = Hash.new
		enabled = false
		["o0", "o1", "o2", "a0", "a1", "a2", "b0", "b1", "b2", "c0", "c1", "c2"].each { |k|
		  begin
		    s = value(k)
			if s == nil || s == ""
			  h[k] = 0.0
			else
			  enabled = true
		      h[k] = Float(eval(s))
			  set_value(k, h[k].to_s)
			end
		  rescue
		    mes = "Cannot evaluate #{value(k)}: " + $!.to_s
			Dialog.run("Value Error") {
			  layout(1, item(:text, :title=>mes))
			  set_attr(1, :hidden=>true)
			}
			return nil
		  end
		}
		if enabled
		  ax = Vector3D[h["a0"], h["a1"], h["a2"]]
		  bx = Vector3D[h["b0"], h["b1"], h["b2"]]
		  cx = Vector3D[h["c0"], h["c1"], h["c2"]]
		  ox = Vector3D[h["o0"], h["o1"], h["o2"]]
		  fx = [1, 1, 1]
		  if ax.length2 < 1e-8
		    fx[0] = 0
		    ax = Vector3D[1, 0, 0]
		  end
		  if bx.length2 < 1e-8
		    fx[1] = 0
		    bx = Vector3D[0, 1, 0]
		  end
		  if cx.length2 < 1e-8
		    fx[2] = 0
		    cx = Vector3D[0, 0, 1]
		  end
		  @box = [ax, bx, cx, ox, fx]
		else
		  @box = nil
		end
		update_values(item1)
	  end
	  def clear_box(item1)
	    @box_save = @box
		@box = nil
		set_attr("restore", :enabled=>true)
		update_values(item1)
	  end
	  def restore_box(item1)
	    @box = @box_save.dup
		set_attr("restore", :enabled=>false)
		update_values(item1)
	  end
	  def angle(v1, v2)
	    acos(v1.dot(v2)/(v1.length * v2.length)) * 180.0 / PI
	  end
	  def update_values(it)
	    ["a0", "a1", "a2", "b0", "b1", "b2", "c0", "c1", "c2", "o0", "o1", "o2",
		 "a", "b", "c", "alpha", "beta", "gamma"].each_with_index { |tag, i|
		  if @box == nil
		    set_value(tag, "")
		  else
		    if i < 12
		      val = @box[i / 3][i % 3]
			  fmt = "%.8g"
		    elsif i < 15
		      val = @box[i - 12].length
			  fmt = "%g"
		    else
		      val = angle(@box[(i - 14) % 3], @box[(i - 13) % 3])
			  fmt = "%g"
			end
			set_value(tag, sprintf(fmt, val))
		  end
		}
	  end

	  layout(4,
	    item(:text, :title=>"Cell parameters:"),
		-1, -1, -1,
		
		item(:text, :title=>"a/b/c"),
		item(:textfield, :width=>140, :tag=>"a"),
		item(:textfield, :width=>140, :tag=>"b"),
		item(:textfield, :width=>140, :tag=>"c"),
		
		item(:text, :title=>"α/β/γ"),
		item(:textfield, :width=>140, :tag=>"alpha"),
		item(:textfield, :width=>140, :tag=>"beta"),
		item(:textfield, :width=>140, :tag=>"gamma"),

		-1, -1, -1,
        item(:button, :title=>"Set Cell Parameters", :align=>:right, :action=>:set_cell_value),
        item(:line),
		-1, -1, -1,

	    item(:text, :title=>"Axis vectors:"),
		-1, -1, -1,

	    item(:text, :title=>"origin"),
		item(:textfield, :width=>140, :tag=>"o0"),
		item(:textfield, :width=>140, :tag=>"o1"),
		item(:textfield, :width=>140, :tag=>"o2"),

	    item(:text, :title=>"a-axis"),
		item(:textfield, :width=>140, :tag=>"a0"),
		item(:textfield, :width=>140, :tag=>"a1"),
		item(:textfield, :width=>140, :tag=>"a2"),

	    item(:text, :title=>"b-axis"),
		item(:textfield, :width=>140, :tag=>"b0"),
		item(:textfield, :width=>140, :tag=>"b1"),
		item(:textfield, :width=>140, :tag=>"b2"),

	    item(:text, :title=>"c-axis"),
		item(:textfield, :width=>140, :tag=>"c0"),
		item(:textfield, :width=>140, :tag=>"c1"),
		item(:textfield, :width=>140, :tag=>"c2"),

        -1, -1, -1,
        item(:button, :title=>"Set Axis Vectors", :align=>:right, :action=>:set_box_value),
        item(:line),
		-1, -1, -1,
	
		item(:text, :title=>"(Ruby expressions are allowed as the values; e.g. 12.0*cos(60*PI/180))"),
		-1, -1, -1,

		item(:button, :title=>"Clear Cell", :tag=>"clear", :action=>:clear_box),
		item(:button, :title=>"Restore Cell", :tag=>"restore", :action=>:restore_box, :enabled=>false),
		-1, -1)
	  set_attr(0, :action=>lambda { |it|
	    if set_box_value(it)
		  new_box = @box
		  end_modal(it)
		end
	  } )
	  update_values(nil)
	}
	if hash[:status] == 0
	  if new_box != nil
	    h = Dialog.run(" ") {
		  layout(1,
		    item(:text, :title=>"Do you want to keep the fractional coordinates?"),
			item(:radio, :title=>"Yes, convert the cartesian coordinates so that the fractional coordinates will be unchanged.",
			  :tag=>"yes", :value=>1),
			item(:radio, :title=>"No, keep the cartesian coordinates unchanged.", :tag=>"no"))
		  radio_group("yes", "no")
		}
		if h[:status] == 0
		  yes = h["yes"] != 0
		  new_box.push(yes)
		  mol.box = new_box
		end
	  else
	    mol.box = nil
	  end
	end
  end


#  For debug; check the symmetry data in space_groups.txt for self-consistency
def Molecule.check_space_group
  zero_transform = Transform.zero
  group_member_p = lambda { |g, tr|
	g.each_with_index { |tr2, idx|
	  tr2 = tr2 - tr
	  #  Wrap the translational part to [-0.5..0.5]
	  3.times { |k| tr2[3, k] -= (tr2[3, k] + 0.5).floor }
	  if tr2.nearly_equal?(zero_transform)
		return idx
	  end
	}
	return nil
  }
  @@space_groups.each { |g|
    name = g[0]
	group = g[1]
	n = group.count
	err = ""
	n.times { |i|
	  tri = group[i]
	  begin
	    trx = tri.inverse
	  rescue
	    err += "The element #{i + 1} is not regular transform; #{trx ? symmetry_to_string(trx) : 'nil'}\n"
		next
	  end
	  3.times { |k| trx[3, k] -= trx[3, k].floor }  #  Wrap the translation part to [0..1]
	  if !group_member_p.call(group, trx)
		err += "The element #{i + 1} has no inverse element; #{symmetry_to_string(trx)}\n"
		next
	  end
	  n.times { |j|
		trj = group[j]
		trx = tri * trj
		3.times { |k| trx[3, k] -= trx[3, k].floor }  #  Wrap the translation part to [0..1]
		if !group_member_p.call(group, trx)
		  err += "The element #{i + 1} * #{j + 1} is not in the group; "
		  err += "#{symmetry_to_string(tri)} * #{symmetry_to_string(trj)} = #{symmetry_to_string(trx)}\n"
		end
	  }
	}
	if err == ""
	  msg = "#{name} is a complete group\n"
	else
	  msg = "#{name} is NOT a complete group\n" + err
	  error_message_box(msg)
	end
#	break if err != ""
  }
end

def Molecule.setup_space_groups
  if defined?(@@space_groups)
    return @@space_groups
  end

  #  The debug flag causes self-check of the space group table
  debug = false
  
  lines = [] if debug

  @@space_groups = []       #  Sequential table of space groups (array of [name, array_of_transform])
  @@space_group_list = []   #  [["Triclinic",    -> Crystal system
							#     ["#1: P1",     -> Space group family
							#       [["P1", 0]]],-> Variations
							#     ...],
							#   ["Monoclinic",
							#     ...
							#     ["#7: Pc",
							#       [["Pc", 6], ["Pa", 7], ["Pn", 8]]],
							#     ...]]
  last_group_code = 0
  Kernel.open(ScriptPath + "/space_groups.txt", "r") { |fp|
	while ln = fp.gets
	  lines.push(ln) if debug
	  if ln =~ /\[(\w+)\]/
		label = $1
		@@space_group_list.push([label, []])
		next
	  end
	  if ln =~ /^(\d+) +(.*)$/
		group_code = $1.to_i
		group_name = $2
	    name = "#" + $1 + ": " + $2
		group = []
		group_gen = []
		while ln = fp.gets
		  lines.push(ln.dup) if debug
		  ln.chomp!
		  if ln =~ /^\s*$/
			break
		  elsif ln =~ /^\+\(/
			lattice = ln.scan(/((^\+)|,)\(([\/,0-9]+)\)/).map {
			  |a| a[2].split(/, */).map {
				|x| if x =~ /^(\d+)\/(\d+)$/
				      $1.to_f / $2.to_f
					else
					  x.to_f
					end } }
			lattice.each { |lat|
			  group.each { |tr|
				tr = tr.dup
				tr[3, 0] += lat[0]
				tr[3, 1] += lat[1]
				tr[3, 2] += lat[2]
				group_gen.push(tr)
			  }
			}
		  else
			#  Rectify the shorthand descriptions
			ln1 = ln.dup if debug
			mod = false
			mod = ln.gsub!(/(\d)(\d)/, "\\1/\\2")
			mod = ln.gsub!(/(^|[^\/])(\d)([^\/]|$)/, "\\11/\\2\\3") || mod
			mod = ln.gsub!(/([0-9xyzXYZ])([xyzXYZ])/, "\\1+\\2") || mod
			if debug && mod
			  puts "#{fp.lineno}: #{ln1} -> #{ln}\n"
			end
			begin
			  group.push(string_to_symmetry(ln))
			  if debug
				lines[-1] = symmetry_to_string(group[-1]).gsub(/ +/, "") + "\n"
			  end
			rescue
			  error_message_box($!.message + "(line #{fp.lineno})")
			end
		  end
		end
		group.concat(group_gen)
		@@space_groups.push([name, group])
		group_index = @@space_groups.count - 1
		if last_group_code != group_code
		  idx = (name.index(" (") || 0) - 1  #  Remove parenthesized phrase
		  group_family_name = name[0..idx]
		  @@space_group_list[-1][1].push([group_family_name, []])
		  last_group_code = group_code
		end
		@@space_group_list[-1][1][-1][1].push([group_name, group_index])
	  end
	end
  }
  if debug
	Molecule.check_space_group
	fn = Dialog.save_panel("Save modified space_group.txt:")
	if fn
	  Kernel.open(fn, "w") { |fp|
		lines.each { |ln|
		  fp.print ln
		}
	  }
	end
  end
  return @@space_groups
end

def cmd_symmetry_operations

  mol = self
  syms = mol.symmetries

  if !defined?(@@space_groups)
    Molecule.setup_space_groups
  end

  #  Menu titles
  crystal_systems = @@space_group_list.map { |m| m[0] }
  space_groups = @@space_group_list[0][1].map { |m| m[0] }
  variations = @@space_group_list[0][1][0][1].map { |m| m[0] }
  
  #  The index (for @@space_groups) of the current space group (an integer or nil)
  current_space_group = nil

  #  Find current space group
  guess_space_group = lambda {
    current_space_group = nil
    s = syms.map { |tr| tr = tr.dup; 3.times { |k| tr[3, k] -= tr[3, k].floor }; tr }
	@@space_groups.each_with_index { |g, i|
	  next if g[1].count != s.count
	  ss = s.dup
	  g[1].each { |tr|
	    idx = ss.find_index { |tr2| tr.nearly_equal?(tr2) }
		break if idx == nil
		ss.delete_at(idx)
	  }
	  if ss.empty?
	    current_space_group = i
		break
	  end
	}
  }
  
  select_space_group = lambda { |it|

    h = Dialog.run("Select Space Group") {
	  update_operation = lambda { |*d|
	    it = item_with_tag("operations")
		cval = item_with_tag("crystal_system")[:value]
		sval = item_with_tag("space_group")[:value]
		vval = item_with_tag("variation")[:value]
		idx = @@space_group_list[cval][1][sval][1][vval][1]
	    op = ""
	    @@space_groups[idx][1].each { |tr|
	      op += symmetry_to_string(tr) + "\n"
	    }
	    it[:value] = op
	  }
      crystal_system_popup_handler = lambda { |it|
	    val = it[:value]
        space_groups = @@space_group_list[val][1].map { |m| m[0] }
        variations = @@space_group_list[val][1][0][1].map { |m| m[0] }
		item_with_tag("space_group")[:subitems] = space_groups
		item_with_tag("space_group")[:value] = 0
		item_with_tag("variation")[:subitems] = variations
		item_with_tag("variation")[:value] = 0
		update_operation.call
      }
      space_group_popup_handler = lambda { |it|
	    val = it[:value]
		cval = item_with_tag("crystal_system")[:value]
		variations = @@space_group_list[cval][1][val][1].map { |m| m[0] }
		item_with_tag("variation")[:subitems] = variations
		item_with_tag("variation")[:value] = 0
		update_operation.call
      }
	  variation_popup_handler = lambda { |it|
	    update_operation.call
	  }
	  layout(2,
	    item(:text, :title=>"Crystal System"),
		item(:popup, :subitems=>crystal_systems, :tag=>"crystal_system", :action=>crystal_system_popup_handler),
		item(:text, :title=>"Space Group"),
		item(:popup, :subitems=>space_groups, :tag=>"space_group", :action=>space_group_popup_handler),
		item(:text, :title=>"Variation"),
		item(:popup, :subitems=>variations, :tag=>"variation", :action=>variation_popup_handler),
		item(:textview, :width=>360, :height=>200, :editable=>false, :tag=>"operations"),
		-1)
	  if current_space_group
	    catch(:end_loop) {
	      @@space_group_list.each_with_index { |m, i|
		    m[1].each_with_index { |mm, j|
		      mm[1].each_with_index { |mmm, k|
		  	    if mmm[1] == current_space_group
			      (it = item_with_tag("crystal_system"))[:value] = i
				  crystal_system_popup_handler.call(it)
			      (it = item_with_tag("space_group"))[:value] = j
				  space_group_popup_handler.call(it)
			      (it = item_with_tag("variation"))[:value] = k
				  variation_popup_handler.call(it)
				  throw(:end_loop)
			    end
			  }
		    }
		  }
		}
	  end
	}
	if h[:status] == 0
	  cval = h["crystal_system"]
	  sval = h["space_group"]
	  vval = h["variation"]
	  idx = @@space_group_list[cval][1][sval][1][vval][1]
	  syms = @@space_groups[idx][1].dup
	end
  }
  
  h = Dialog.run("Symmetry Operations") {
    update_space_group_text = lambda {
	  guess_space_group.call
	  if current_space_group
	    label = @@space_groups[current_space_group][0]
	  else
	    label = ""
	  end
	  item_with_tag("space_group")[:value] = label
	}
    layout(1,
	  layout(2,
	    item(:text, :title=>"Space Group:"),
		item(:textfield, :width=>200, :editable=>false, :tag=>"space_group")),
	  item(:button, :title=>"Select...", :action=>lambda { |it|
		select_space_group.call(it)
		item_with_tag("sym_table")[:refresh] = true
		update_space_group_text.call }, :align=>:right),
	  item(:table, :tag=>"sym_table", :width=>300, :height=>300,
		:columns=>[["Symmetry Operations", 280]],
        :on_count=>lambda { |it| syms.count },
        :on_get_value=>lambda { |it, row, col| symmetry_to_string(syms[row]) },
        :on_set_value=>lambda { |it, row, col, val|
		  syms[row] = string_to_symmetry(val)
		  update_space_group_text.call },
        :is_item_editable=>lambda { |it, row, col| true },
        :on_selection_changed=>lambda { |it|
		  item_with_tag("delete")[:enabled] = (it[:selection].count > 0) } ),
	  layout(2,
	    item(:togglebutton, :title=>"+", :width=>40, :tag=>"add",
		  :action=>lambda { |it|
		    syms.push(Transform.new); item_with_tag("sym_table")[:refresh] = true
			update_space_group_text.call } ),
		item(:togglebutton, :title=>"-", :width=>40, :tag=>"delete", :enabled=>false,
		  :action=>lambda { |it|
		    item_with_tag("sym_table")[:selection].reverse_each { |n| syms.delete_at(n) }
			item_with_tag("sym_table")[:refresh] = true
			update_space_group_text.call } ),
		:padding=>0, :margin=>0 ) )
	item_with_tag("sym_table")[:refresh] = true
	update_space_group_text.call
  }
  if h[:status] == 0
    self.remove_symmetries(nil)
	syms.each { |s|
	  self.add_symmetry(s)
    }
  end
end

def cmd_show_periodic_image
    mol = self
    hash = Dialog.run("Show Periodic Image") {
	  @mol = mol
	  def set_periodic_image(it)
	    a = []
	    ["amin", "amax", "bmin", "bmax", "cmin", "cmax"].each_with_index { |k, i|
		  s = value(k)
		  if s == nil || s == ""
		    a[i] = 0
		  else
		    a[i] = Integer(s)
		  end
		}
	    @mol.show_periodic_image(a)
		@mol.show_periodic_image = (attr("show_flag", :value) == 1 ? true : false)
	  end
	  pimage = @mol.show_periodic_image
	  flag = @mol.show_periodic_image?
	  layout(4,
	    item(:checkbox, :title=>"Show Periodic Image", :tag=>"show_flag", :value=>(flag ? 1 : 0),
		  :action=>lambda { |it| @mol.show_periodic_image = (it[:value] == 1 ? true : false) } ),
		-1, -1, -1,
		item(:text, :title=>"a-axis"),
		item(:textfield, :width=>80, :tag=>"amin", :value=>pimage[0].to_s),
		item(:text, :title=>"to"),
		item(:textfield, :width=>80, :tag=>"amax", :value=>pimage[1].to_s),
		item(:text, :title=>"b-axis"),
		item(:textfield, :width=>80, :tag=>"bmin", :value=>pimage[2].to_s),
		item(:text, :title=>"to"),
		item(:textfield, :width=>80, :tag=>"bmax", :value=>pimage[3].to_s),
		item(:text, :title=>"c-axis"),
		item(:textfield, :width=>80, :tag=>"cmin", :value=>pimage[4].to_s),
		item(:text, :title=>"to"),
		item(:textfield, :width=>80, :tag=>"cmax", :value=>pimage[5].to_s),
		item(:button, :title=>"Set", :action=>:set_periodic_image))
	  set_attr(0, :action=>lambda { |it| set_periodic_image(it); end_modal(it) } )
	}
end

def has_expanded_atoms
  atoms.find { |ap| ap.symop != nil }
end

def cmd_remove_expanded_atoms
  mol = self
  return unless mol.has_expanded_atoms
  hash = Dialog.run("Remove Expanded Atoms") {
    layout(1,
	  item(:radio, :title=>"Remove all expanded atoms.", :tag=>"all", :value=>1),
	  item(:radio, :title=>"Keep the expanded atoms that are in the same fragment with original atoms.", :tag=>"partial"))
	radio_group("all", "partial")
  }
  if hash[:status] == 0
    if hash["all"] == 1
	  rm = mol.atom_group { |ap| ap.symop != nil }
	else
	  rm = IntGroup[]
	  mol.each_fragment { |f|
	    if !f.find { |i| mol.atoms[i].symop == nil }
		  rm.add(f)
		end
	  }
	end
	mol.remove(rm)
  end
end

end

require_cell = lambda { |m| m && m.cell != nil }

register_menu("Xtal\tUnit Cell...", :cmd_unit_cell)
register_menu("Xtal\tSymmetry Operations...", :cmd_symmetry_operations)
register_menu("Xtal\t-", nil)
register_menu("Xtal\tComplete by Symmetry", :complete_by_symmetry, lambda { |m| m && (m.cell != nil || m.nsymmetries > 1) } )
register_menu("Xtal\tCreate Packing Diagram...", :create_packing_diagram, require_cell)
register_menu("Xtal\tShow Periodic Image...", :cmd_show_periodic_image, require_cell)
register_menu("Xtal\tRemove Expanded Atoms...", :cmd_remove_expanded_atoms, lambda { |m| m && m.has_expanded_atoms } )
register_menu("Xtal\t-", nil)
register_menu("Xtal\tBonds and Angles with Sigma...", :cmd_bond_angle_with_sigma, :non_empty)
register_menu("Xtal\tBest-fit Planes...", :cmd_plane, :non_empty)
register_menu("Xtal\tShow ORTEP...", :cmd_show_ortep, :non_empty)
