#!/usr/bin/env python

from lxml import etree
from pylab import *

ns = {"cml":"http://www.xml-cml.org/schema"}

doc = etree.parse("siesta.xml")

molecule = doc.xpath("//cml:molecule[1]", ns)[0]

elementNames = molecule.xpath("cml:atomArray/cml:atom/@elementType", ns)

print elementNames

numAtoms = molecule.xpath("count(cml:atomArray/cml:atom)", ns)
numCls = molecule.xpath("count(cml:atomArray/cml:atom[@elementType='Cl'])", ns)

print 'This molecule has ', numCls, ' chlorines.'

xCoords = molecule.xpath("cml:atomArray/cml:atom/@x3", ns)
yCoords = molecule.xpath("cml:atomArray/cml:atom/@y3", ns)
zCoords = molecule.xpath("cml:atomArray/cml:atom/@z3", ns)

# These next three lines convert the value retrieved from the
# XML file (which is always a string) to a number for
# manipulation below. This is a Python way of doing it - some
# other languages may not need this.
xCoords = [float(x) for x in xCoords]
yCoords = [float(y) for y in yCoords]
zCoords = [float(z) for z in zCoords]

xCentre = 0
yCentre = 0
zCentre = 0
for i in range(int(numAtoms)):
    xCentre += xCoords[i]
    yCentre += yCoords[i]
    zCentre += zCoords[i]

xCentre /= numAtoms
yCentre /= numAtoms
zCentre /= numAtoms

print 'The geometric centre of the molecule is at ', xCentre, yCentre, zCentre


