#!@PYTHON@
#
# CreateVspeShapes  v1.0 by neuen 07May 2009
#
# This programm creates basic VSEPR Structures
# Given an element as input this programm obtains the Valence number of the Atom and
# creates a molecule with that element and as many H atoms bonded to it as possible.
# Between 0 and 12 it will give the molecule CLOSE to optimal shape, optimization is HIGHLY advised.
#

import sys, math
werr = sys.stderr.write

Bondlength = 1.0
temp_Bondlength = 0.0
ValenceNumber = 0
vector = [0.0, 0.0, 0.0]
if len(sys.argv)<3:
  werr("Usage is " + sys.argv[0] + " <element abbreviation> <required valence number>\n")
  sys.exit(1)

element = sys.argv[1]
ValenceNumber = int(sys.argv[2])

# get valence number from file

output = open(element+"H_%d.xyz" %(ValenceNumber), "w")

if (ValenceNumber > 12):
  werr("the Valence Number of this element is larger than 12, the corresponding shape has not been implemented! \n")
elif (ValenceNumber < 0):
  werr("The Valence Number of this element is smaller than 0, something is wrong with the database! \n")

else:
#Write out header of xyz file and single atom of element
  output.write( str(ValenceNumber+1) +" \n \n")
  output.write(element + " 0.0 0.0 0.0 \n")
  if (ValenceNumber >= 1):
    # We will always need another molecule anyplace, all the others orient themselves depending on it.
    vector = [Bondlength, 0.0, 0.0]
    output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
    
    if (ValenceNumber == 2 or (ValenceNumber >4 and ValenceNumber < 8)):
      # In those cases we always get a line & a circle with the rest.
      vector = [-Bondlength, 0.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      if (ValenceNumber == 5):
        vector = [0.0, Bondlength, 0.0]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, -Bondlength*0.5, Bondlength*math.sqrt(3.0)*0.5]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, -Bondlength*0.5, -Bondlength*0.5*math.sqrt(3.0)]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      
      elif (ValenceNumber == 6):
        vector = [0.0, Bondlength, 0.0]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, -Bondlength, 0.0]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, 0.0, Bondlength]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, 0.0, -Bondlength]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        
      elif (ValenceNumber == 7):
        vector = [0.0, Bondlength, 0.0]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, Bondlength*math.cos(math.pi*0.4), Bondlength*math.sin(math.pi*0.4)]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, Bondlength*math.cos(math.pi*0.8), Bondlength*math.sin(math.pi*0.8)]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, Bondlength*math.cos(math.pi*1.2), Bondlength*math.sin(math.pi*1.2)]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        vector = [0.0, Bondlength*math.cos(math.pi*1.6), Bondlength*math.sin(math.pi*1.6)]
        output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
        
    elif (ValenceNumber == 3):
      vector = [-Bondlength*0.5, Bondlength*math.sqrt(3)*0.5, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength*0.5, -Bondlength*math.sqrt(3)*0.5, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

    elif (ValenceNumber == 4):
      vector = [-Bondlength/3.0, Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

    elif (ValenceNumber == 8):
      vector = [-Bondlength, 0.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength/3.0, Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

      vector = [Bondlength/3.0, -Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength/3.0, Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength/3.0, Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

    

    elif (ValenceNumber ==9):
      vector = [Bondlength*math.cos(0.4*math.pi), Bondlength*math.sin(0.4*math.pi), 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), -Bondlength*math.sin(0.4*math.pi), 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), 0.0, Bondlength*math.sin(0.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), 0.0, -Bondlength*math.sin(0.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

      vector = [Bondlength*math.cos(0.8*math.pi), Bondlength*math.sin(0.8*math.pi)*math.sin(0.25*math.pi), Bondlength*math.sin(0.8*math.pi)*math.cos(0.25*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), Bondlength*math.sin(0.8*math.pi)*math.sin(0.75*math.pi), Bondlength*math.sin(0.8*math.pi)*math.cos(0.75*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), Bondlength*math.sin(0.8*math.pi)*math.sin(1.25*math.pi), Bondlength*math.sin(0.8*math.pi)*math.cos(1.25*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), Bondlength*math.sin(0.8*math.pi)*math.sin(1.75*math.pi), Bondlength*math.sin(0.8*math.pi)*math.cos(1.75*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      

    elif (ValenceNumber ==10):
      vector = [-Bondlength, 0.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)

      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), 0.0, temp_Bondlength]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), 0.0, -temp_Bondlength]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(0.25*math.pi), temp_Bondlength*math.cos(0.25*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(0.75*math.pi), temp_Bondlength*math.cos(0.75*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(1.25*math.pi), temp_Bondlength*math.cos(1.25*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(1.75*math.pi), temp_Bondlength*math.cos(1.75*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      


    elif (ValenceNumber == 11):
      temp_Bondlength = Bondlength*math.sin(0.4*math.pi)
      #above tempral value will be used in order to save sin calculations in every atom
      vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(0.4*math.pi),  temp_Bondlength*math.sin(0.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(0.8*math.pi),  temp_Bondlength*math.sin(0.8*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(1.2*math.pi),  temp_Bondlength*math.sin(1.2*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(1.6*math.pi),  temp_Bondlength*math.sin(1.6*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      temp_Bondlength = Bondlength*math.sin(0.8*math.pi)
      #above tempral value will be used in order to save sin calculations in every atom
      vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(0.2*math.pi),  temp_Bondlength*math.sin(0.2*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(0.6*math.pi),  temp_Bondlength*math.sin(0.6*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.0*math.pi),  temp_Bondlength*math.sin(1.0*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.4*math.pi),  temp_Bondlength*math.sin(1.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.8*math.pi),  temp_Bondlength*math.sin(1.8*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      
      

    elif (ValenceNumber == 12):
      vector = [-Bondlength, 0.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)
      # above temporal value will save a lot of divisions
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(0.4*math.pi),  temp_Bondlength*math.sin(0.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(0.8*math.pi),  temp_Bondlength*math.sin(0.8*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(1.2*math.pi),  temp_Bondlength*math.sin(1.2*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(1.6*math.pi),  temp_Bondlength*math.sin(1.6*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

      vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(0.4*math.pi),  -temp_Bondlength*math.sin(0.4*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(0.8*math.pi),  -temp_Bondlength*math.sin(0.8*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(1.2*math.pi),  -temp_Bondlength*math.sin(1.2*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(1.6*math.pi),  -temp_Bondlength*math.sin(1.6*math.pi)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      
# Cases below will never be called as the program is (all inquiries larger 12 blocked)
# However I wrote 14 erronously instead of 12, I leave this here, should it be needed at some point.

    elif (ValenceNumber == 13):
      werr("This ValenceNumber has not been implemented\n.")

    elif (ValenceNumber == 14):
      vector = [-Bondlength, 0.0, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)
      # above temporal value will save a lot of divisions
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(math.pi/3.0), temp_Bondlength*math.sin(math.pi/3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(2.0*math.pi/3.0), temp_Bondlength*math.sin(2.0*math.pi/3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(4.0*math.pi/3.0), temp_Bondlength*math.sin(4.0*math.pi/3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(5.0*math.pi/3.0), temp_Bondlength*math.sin(5.0*math.pi/3.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")

      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(math.pi/6.0),  temp_Bondlength*math.sin(math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(math.pi/3.0 +math.pi/6.0),  temp_Bondlength*math.sin(math.pi/3.0+math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(2.0*math.pi/3.0 +math.pi/6.0),  temp_Bondlength*math.sin(2.0*math.pi/3.0 +math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(3.0*math.pi/3.0 +math.pi/6.0),  temp_Bondlength*math.sin(3.0*math.pi/3.0 +math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(4.0*math.pi/3.0 +math.pi/6.0),  temp_Bondlength*math.sin(4.0*math.pi/3.0 +math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(5.0*math.pi/3.0 +math.pi/6.0),  temp_Bondlength*math.sin(5.0*math.pi/3.0 +math.pi/6.0)]
      output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
      

output.close()