#!/usr/bin/python

import System
from math import sqrt


# a and c are lattice parameters for hex. ## I think.  double check "s3"
s3=sqrt(3)
a= s3
c= (5.27/3.31)*a
m=7.4651033e-26

cell=[
  a*sqrt(3)/2,-a/2,0,
  0,a,0,
  0,0,c,
]

atoms=[
  [ "A", m ],
  [ "B", m ],
]

sites=[
  [ 0,0,0,             0 ],
  [ a/s3,0,c/2,     1 ],
]

bonds=[
  [ 0,1, System.modAxial([ a/s3,0, c/2],   23.198,[-2.576,-2.576,-2.559]) ],
  [ 1,0, System.modAxial([-a/s3,0,-c/2],   23.198,[-2.576,-2.576,-2.559]) ],

  [ 0,0, System.modAxial([ 0, a,0],          11.793,[ 2.374, 2.374, 1.525]) ],
  [ 1,1, System.modAxial([ 0, a,0],          11.793,[ 2.374, 2.374, 1.525]) ],

  [ 0,1, System.modAxial([-2*a/s3,0, c/2], -9.944,[ 2.284, 2.284, 1.981]) ],
  [ 1,0, System.modAxial([ 2*a/s3,0,-c/2], -9.944,[ 2.284, 2.284, 1.981]) ],

  [ 0,0, System.modAxial([ 0,0,c],            0.00,[-0.622,-0.622, -0.610]) ],
  [ 1,1, System.modAxial([ 0,0,c],            0.00,[-0.622,-0.622, -0.610]) ],

  [ 0,1, System.modAxial([ 5*a/(2*s3), a/2, c/2],  
                                                3.557,[ 0.539, 0.539,-0.648]) ],
  [ 1,0, System.modAxial([-5*a/(2*s3),-a/2,-c/2],  
                                                3.557,[ 0.539, 0.539,-0.648]) ],

  [ 0,0, System.modAxial([ a*s3,0,0],        1.307,[-0.041,-0.041,-0.115]) ],
  [ 1,1, System.modAxial([ a*s3,0,0],        1.307,[-0.041,-0.041,-0.115]) ],
  [ 0,0, System.modAxial([-a*s3,0,0],        1.307,[-0.041,-0.041,-0.115]) ],
  [ 1,1, System.modAxial([-a*s3,0,0],        1.307,[-0.041,-0.041,-0.115]) ],

]

System.write(cell,atoms,sites,bonds,"hcp")
