[MMTK] Re: Defining ciclobutadiene

[Konrad Hinsen]

> I'm sorry to be such a pain in the rear but I have corrected the
> Z-Matrix:
> 
...

> and I get a *very* similar error message:
>
...
>   File "/usr/local/mmtk/ConfigIO.py", line 169, in findPositions
>     pos3 = self.coordinates[entry[5]]
> KeyError: <Atom number 3>

This time you found a bug in MMTK (and in one of its oldest parts!).
The geometrical construction of the atom positions doesn't work
for your rather special case with plenty of right angles and
zero dihedrals.

To fix the bug, replace the class ZMatrix in the file ConfigIO.py
by the following one:

class ZMatrix:

    def __init__(self, data):
	self.data = data
	self.coordinates = {}

    substitute = 1

    def findPositions(self):
	# First atom at origin
	self.coordinates[self.data[0][0]] = Vector(0,0,0)
	# Second atom along x-axis
	self.coordinates[self.data[1][0]] = Vector(self.data[1][2],0,0)
	# Third atom in xy-plane
        try:
            pos1 = self.coordinates[self.data[2][1]]
        except KeyError:
            raise ValueError, \
                  "atom %d has no defined position" % self.data[2][1].number
        try:
            pos2 = self.coordinates[self.data[2][3]]
        except KeyError:
            raise ValueError, \
                  "atom %d has no defined position" % self.data[2][3].number
	sphere = Sphere(pos1, self.data[2][2])
	cone = Cone(pos1, pos2-pos1, self.data[2][4])
	plane = Plane(Vector(0,0,0), Vector(0,0,1))
	points = sphere.intersectWith(cone).intersectWith(plane)
	self.coordinates[self.data[2][0]] = points[0]
	# All following atoms defined by distance + angle + dihedral
	for entry in self.data[3:]:
            try:
                pos1 = self.coordinates[entry[1]]
            except KeyError:
                raise ValueError, \
                      "atom %d has no defined position" % entry[1].number
            try:
                pos2 = self.coordinates[entry[3]]
            except KeyError:
                raise ValueError, \
                      "atom %d has no defined position" % entry[3].number
            try:
                pos3 = self.coordinates[entry[5]]
            except KeyError:
                raise ValueError, \
                      "atom %d has no defined position" % entry[5].number
	    distance = entry[2]
	    angle = entry[4]
	    dihedral = entry[6]
	    sphere = Sphere(pos1, distance)
	    cone = Cone(pos1, pos2-pos1, angle)
	    plane123 = Plane(pos3, pos2, pos1)
            points = sphere.intersectWith(cone).intersectWith(plane123)
            for p in points:
                if Plane(pos2, pos1, p).normal * plane123.normal > 0:
                    break
            p = rotatePoint(p, Line(pos1, pos2-pos1), dihedral)
            self.coordinates[entry[0]] = p

    def applyTo(self, object):
	if not len(self.coordinates):
	    self.findPositions()
	for entry in self.coordinates.items():
	    object.setPosition(entry[0], entry[1])
	object.normalizePosition()

With that definition I can construct your molecule; I think you still
need to correct one angle to make it come out right.

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hinsen@cnrs-orleans.fr
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.55.69
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------