source: util/src/CreateVspeShapes.py.in@ aa8542

Last change on this file since aa8542 was e1a46d, checked in by Frederik Heber <heber@…>, 16 years ago

All of my python script put into ESPACK and adapted with @PYTHON@ and so on.
  • Property mode set to 100755
File size: 17.2 KB
RevLine 
[e1a46d]1#!@PYTHON@
2#
3# CreateVspeShapes v1.0 by neuen 07May 2009
4#
5# This programm creates basic VSEPR Structures
6# Given an element as input this programm obtains the Valence number of the Atom and
7# creates a molecule with that element and as many H atoms bonded to it as possible.
8# Between 0 and 12 it will give the molecule CLOSE to optimal shape, optimization is HIGHLY advised.
9#
10
11import sys, math
12werr = sys.stderr.write
13
14Bondlength = 1.0
15temp_Bondlength = 0.0
16ValenceNumber = 0
17vector = [0.0, 0.0, 0.0]
18if len(sys.argv)<3:
19 werr("Usage is " + sys.argv[0] + " <element abbreviation> <required valence number>\n")
20 sys.exit(1)
21
22element = sys.argv[1]
23ValenceNumber = int(sys.argv[2])
24
25# get valence number from file
26
27output = open(element+"H_%d.xyz" %(ValenceNumber), "w")
28
29if (ValenceNumber > 12):
30 werr("the Valence Number of this element is larger than 12, the corresponding shape has not been implemented! \n")
31elif (ValenceNumber < 0):
32 werr("The Valence Number of this element is smaller than 0, something is wrong with the database! \n")
33
34else:
35#Write out header of xyz file and single atom of element
36 output.write( str(ValenceNumber+1) +" \n \n")
37 output.write(element + " 0.0 0.0 0.0 \n")
38 if (ValenceNumber >= 1):
39 # We will always need another molecule anyplace, all the others orient themselves depending on it.
40 vector = [Bondlength, 0.0, 0.0]
41 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
42
43 if (ValenceNumber == 2 or (ValenceNumber >4 and ValenceNumber < 8)):
44 # In those cases we always get a line & a circle with the rest.
45 vector = [-Bondlength, 0.0, 0.0]
46 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
47 if (ValenceNumber == 5):
48 vector = [0.0, Bondlength, 0.0]
49 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
50 vector = [0.0, -Bondlength*0.5, Bondlength*math.sqrt(3.0)*0.5]
51 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
52 vector = [0.0, -Bondlength*0.5, -Bondlength*0.5*math.sqrt(3.0)]
53 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
54
55 elif (ValenceNumber == 6):
56 vector = [0.0, Bondlength, 0.0]
57 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
58 vector = [0.0, -Bondlength, 0.0]
59 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
60 vector = [0.0, 0.0, Bondlength]
61 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
62 vector = [0.0, 0.0, -Bondlength]
63 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
64
65 elif (ValenceNumber == 7):
66 vector = [0.0, Bondlength, 0.0]
67 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
68 vector = [0.0, Bondlength*math.cos(math.pi*0.4), Bondlength*math.sin(math.pi*0.4)]
69 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
70 vector = [0.0, Bondlength*math.cos(math.pi*0.8), Bondlength*math.sin(math.pi*0.8)]
71 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
72 vector = [0.0, Bondlength*math.cos(math.pi*1.2), Bondlength*math.sin(math.pi*1.2)]
73 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
74 vector = [0.0, Bondlength*math.cos(math.pi*1.6), Bondlength*math.sin(math.pi*1.6)]
75 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
76
77 elif (ValenceNumber == 3):
78 vector = [-Bondlength*0.5, Bondlength*math.sqrt(3)*0.5, 0.0]
79 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
80 vector = [-Bondlength*0.5, -Bondlength*math.sqrt(3)*0.5, 0.0]
81 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
82
83 elif (ValenceNumber == 4):
84 vector = [-Bondlength/3.0, Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
85 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
86 vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
87 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
88 vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
89 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
90
91 elif (ValenceNumber == 8):
92 vector = [-Bondlength, 0.0, 0.0]
93 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
94 vector = [-Bondlength/3.0, Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
95 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
96 vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
97 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
98 vector = [-Bondlength/3.0, -Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
99 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
100
101 vector = [Bondlength/3.0, -Bondlength*2.0*math.sqrt(2.0)/3.0, 0.0]
102 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
103 vector = [Bondlength/3.0, Bondlength*math.sqrt(2.0)/3.0, -Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
104 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
105 vector = [Bondlength/3.0, Bondlength*math.sqrt(2.0)/3.0, Bondlength*math.sqrt(2.0)/math.sqrt(3.0)]
106 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
107
108
109
110 elif (ValenceNumber ==9):
111 vector = [Bondlength*math.cos(0.4*math.pi), Bondlength*math.sin(0.4*math.pi), 0.0]
112 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
113 vector = [Bondlength*math.cos(0.4*math.pi), -Bondlength*math.sin(0.4*math.pi), 0.0]
114 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
115 vector = [Bondlength*math.cos(0.4*math.pi), 0.0, Bondlength*math.sin(0.4*math.pi)]
116 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
117 vector = [Bondlength*math.cos(0.4*math.pi), 0.0, -Bondlength*math.sin(0.4*math.pi)]
118 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
119
120 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)]
121 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
122 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)]
123 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
124 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)]
125 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
126 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)]
127 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
128
129
130 elif (ValenceNumber ==10):
131 vector = [-Bondlength, 0.0, 0.0]
132 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
133 temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)
134
135 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
136 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
137 vector = [temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
138 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
139 vector = [temp_Bondlength/math.sqrt(3.0), 0.0, temp_Bondlength]
140 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
141 vector = [temp_Bondlength/math.sqrt(3.0), 0.0, -temp_Bondlength]
142 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
143
144 vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(0.25*math.pi), temp_Bondlength*math.cos(0.25*math.pi)]
145 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
146 vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(0.75*math.pi), temp_Bondlength*math.cos(0.75*math.pi)]
147 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
148 vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(1.25*math.pi), temp_Bondlength*math.cos(1.25*math.pi)]
149 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
150 vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.sin(1.75*math.pi), temp_Bondlength*math.cos(1.75*math.pi)]
151 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
152
153
154
155 elif (ValenceNumber == 11):
156 temp_Bondlength = Bondlength*math.sin(0.4*math.pi)
157 #above tempral value will be used in order to save sin calculations in every atom
158 vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength, 0.0]
159 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
160 vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.sin(0.4*math.pi)]
161 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
162 vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.sin(0.8*math.pi)]
163 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
164 vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(1.2*math.pi), temp_Bondlength*math.sin(1.2*math.pi)]
165 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
166 vector = [Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.cos(1.6*math.pi), temp_Bondlength*math.sin(1.6*math.pi)]
167 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
168 temp_Bondlength = Bondlength*math.sin(0.8*math.pi)
169 #above tempral value will be used in order to save sin calculations in every atom
170 vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(0.2*math.pi), temp_Bondlength*math.sin(0.2*math.pi)]
171 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
172 vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(0.6*math.pi), temp_Bondlength*math.sin(0.6*math.pi)]
173 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
174 vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.0*math.pi), temp_Bondlength*math.sin(1.0*math.pi)]
175 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
176 vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.4*math.pi), temp_Bondlength*math.sin(1.4*math.pi)]
177 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
178 vector = [Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.cos(1.8*math.pi), temp_Bondlength*math.sin(1.8*math.pi)]
179 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
180
181
182
183 elif (ValenceNumber == 12):
184 vector = [-Bondlength, 0.0, 0.0]
185 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
186 temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)
187 # above temporal value will save a lot of divisions
188 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
189 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
190 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(0.4*math.pi), temp_Bondlength*math.sin(0.4*math.pi)]
191 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
192 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(0.8*math.pi), temp_Bondlength*math.sin(0.8*math.pi)]
193 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
194 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(1.2*math.pi), temp_Bondlength*math.sin(1.2*math.pi)]
195 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
196 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(1.6*math.pi), temp_Bondlength*math.sin(1.6*math.pi)]
197 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
198
199 vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
200 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
201 vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(0.4*math.pi), -temp_Bondlength*math.sin(0.4*math.pi)]
202 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
203 vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(0.8*math.pi), -temp_Bondlength*math.sin(0.8*math.pi)]
204 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
205 vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(1.2*math.pi), -temp_Bondlength*math.sin(1.2*math.pi)]
206 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
207 vector = [-temp_Bondlength/math.sqrt(3.0), -temp_Bondlength*math.cos(1.6*math.pi), -temp_Bondlength*math.sin(1.6*math.pi)]
208 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
209
210# Cases below will never be called as the program is (all inquiries larger 12 blocked)
211# However I wrote 14 erronously instead of 12, I leave this here, should it be needed at some point.
212
213 elif (ValenceNumber == 13):
214 werr("This ValenceNumber has not been implemented\n.")
215
216 elif (ValenceNumber == 14):
217 vector = [-Bondlength, 0.0, 0.0]
218 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
219 temp_Bondlength = Bondlength*0.5*math.sqrt(3.0)
220 # above temporal value will save a lot of divisions
221 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength, 0.0]
222 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
223 vector = [temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(math.pi/3.0), temp_Bondlength*math.sin(math.pi/3.0)]
224 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
225 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)]
226 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
227 vector = [temp_Bondlength/math.sqrt(3.0), -temp_Bondlength, 0.0]
228 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
229 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)]
230 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
231 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)]
232 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
233
234 vector = [-temp_Bondlength/math.sqrt(3.0), temp_Bondlength*math.cos(math.pi/6.0), temp_Bondlength*math.sin(math.pi/6.0)]
235 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
236 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)]
237 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
238 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)]
239 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
240 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)]
241 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
242 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)]
243 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
244 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)]
245 output.write ("H "+ str(vector[0]) +" "+ str(vector[1]) +" "+ str(vector[2]) +" \n")
246
247
248output.close()
Note: See TracBrowser for help on using the repository browser.