import math

class molecule:
    atoms = []
    bonds = []
    def print_atoms(self):
        for atom in self.atoms:
            print(atom)    
    def print_bonds(self):
        """
        >>> moltest.print_bonds()
        1 - 9: 1.228
        2 - 10: 1.236
        3 - 8: 1.388
        3 - 10: 1.417
        3 - 12: 1.458
        4 - 7: 1.365
        4 - 11: 1.371
        4 - 13: 1.443
        5 - 9: 1.420
        5 - 10: 1.442
        5 - 14: 1.459
        6 - 8: 1.369
        6 - 11: 1.317
        7 - 8: 1.368
        7 - 9: 1.427
        11 - 15: 1.082
        12 - 16: 1.094
        12 - 17: 1.093
        12 - 18: 1.094
        13 - 19: 1.093
        13 - 20: 1.093
        13 - 21: 1.093
        14 - 22: 1.092
        14 - 23: 1.095
        14 - 24: 1.095
        """         
        pass
            
    def __init__(self, f):
        self.atoms = []
        self.bonds = []        
        with open(f) as molfile:
            for i in range(4):
                line = molfile.readline()
            na = int(line[0:3])
            nb = 0
            for i in range(na):
                line = molfile.readline()
                x = float(line[0:10])
                y = float(line[10:20])
                z = float(line[20:30])
                symbol = line[31:34].strip()
                self.atoms.append(atom(symbol,x,y,z))
            
class atom:
    symbol = ""
    x = 0.0
    y = 0.0
    z = 0.0
    def __init__(self, symbol, x, y, z):
        self.symbol = symbol
        self.x = x
        self.y = y
        self.z = z
    def __str__(self):
        return("{:3s} {:.3f} {:.3f} {:.3f}".format(self.symbol, self.x, self.y, self.z))    
    def distance(self, a):
        return math.sqrt((self.x-a.x)**2 + (self.y-a.y)**2  + (self.z-a.z)**2)

class bond:
    id1 = 0
    id2 = 0
    bondtype = 0
    def __init__(self, id1, id2, bondtype):
        self.id1 = id1
        self.id2 = id2
        self.bondtype = bondtype
    def distance(self, molecule):
        return molecule.atoms[self.id1-1].distance(molecule.atoms[self.id2-1])

m = molecule("caffeine.mol")
m.print_atoms()
m.print_bonds()

if __name__ == "__main__":
    import doctest
    doctest.testmod(extraglobs={'moltest': molecule("caffeine.mol")})