#!/usr/bin/env python

import math

class Element:
    def __init__(self, name, radius):
        self.name, self.radius = name, radius

elements = {
    "H": Element("Hydrogen", 0.40),
    "C": Element("Carbon",   0.53),
    "N": Element("Nitrogen", 0.53),
    "O": Element("Oxygen",   0.53),
    "F": Element("Fluorine", 0.53),
    }

class Atom:
    def __init__(self, element, pos):
        self.element, self.pos = element, pos

class Molecule:
    def readPdb(self, fp):
        self.name = "Fluoxetine"
        self.atoms = []
        self.bonds = []
        for line in fp.read().split("\r"):
            if line.startswith("ATOM"):
                id, e, dummy, x, y, z = line.split()[1:7]
                assert int(id) == len(self.atoms) + 1
                self.atoms.append(Atom(elements[e], map(float, (x, y, z))))
            elif line.startswith("CONECT"):
                a, b = map(int, line.split()[1:])
                self.bonds.append((self.atoms[a - 1], self.atoms[b - 1]))

    def center(self):
        minpos, maxpos = [0, 0, 0], [0, 0, 0]
        for atom in self.atoms:
            for i in range(3):
                minpos[i] = min(minpos[i], atom.pos[i] - atom.element.radius)
                maxpos[i] = max(maxpos[i], atom.pos[i] - atom.element.radius)
        offsets = [0, 0, 0]
        for i in range(3):
            offsets[i] = (minpos[i] + maxpos[i]) / 2
        for atom in self.atoms:
            for i in range(3):
                atom.pos[i] -= offsets[i]

    def writePov(self, fp):
        bondradius = elements["C"].radius / 2
        print >>fp, '#include "elements.pov"'
        print >>fp, '#declare %s = union {' % self.name
        for atom in self.atoms:
            print >>fp, "  sphere {",
            print >>fp, "<%s, %s, %s>, %s" % (
                tuple(atom.pos) + (atom.element.radius,)),
            print >>fp, "texture { Solid_%s }" % atom.element.name,
            print >>fp, "}"
        for bond in self.bonds:
            if bond[0].element is bond[1].element:
                print >>fp, "  cylinder {"
                print >>fp, "    <%s, %s, %s>," % tuple(bond[0].pos),
                print >>fp, "<%s, %s, %s>," % tuple(bond[1].pos),
                print >>fp, bondradius
                print >>fp, "    texture { Solid_%s }" % bond[0].element.name
                print >>fp, "  }"
            else:
                r = [b - a for a, b in zip(bond[0].pos, bond[1].pos)]
                l = math.sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2])
                k = ((l - bond[0].element.radius - bond[1].element.radius)/2 +
                     bond[0].element.radius) / l
                midpoint = [a + b * k for a, b in zip(bond[0].pos, r)]
                print >>fp, "  cylinder {"
                print >>fp, "    <%s, %s, %s>," % tuple(bond[0].pos),
                print >>fp, "<%s, %s, %s>," % tuple(midpoint),
                print >>fp, bondradius
                print >>fp, "    texture { Solid_%s }" % bond[0].element.name
                print >>fp, "  }"
                print >>fp, "  cylinder {"
                print >>fp, "    <%s, %s, %s>," % tuple(midpoint),
                print >>fp, "<%s, %s, %s>," % tuple(bond[1].pos),
                print >>fp, bondradius
                print >>fp, "    texture { Solid_%s }" % bond[1].element.name
                print >>fp, "  }"
        print >>fp, "  merge {"
        for atom in self.atoms:
            print >>fp, "    sphere {",
            print >>fp, "<%s, %s, %s>, %s" % (
                tuple(atom.pos) + (atom.element.radius * 3,)),
            print >>fp, "texture { Trans_%s }" % atom.element.name,
            print >>fp, "}"
        print >>fp, "  }"
        print >>fp, "}"


if __name__ == "__main__":
    import sys
    mol = Molecule()
    mol.readPdb(sys.stdin)
    mol.center()
    mol.writePov(sys.stdout)