OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
protonSampler
1#!@Python3_EXECUTABLE@
2
3'''
4This python script will (hopefully) generate proton disordered configurations of ice-Ih, given an input of the (.xyz) coordinates of the oxygen positions in the lattice.
5'''
6
7__author__ = "Patrick Louden (plouden@nd.edu)"
8__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
9__license__ = "OpenMD"
10
11import sys
12import argparse
13import textwrap
14import numpy as np
15import random
16import copy
17import math as math
18from fractions import Fraction, gcd
19from argparse import RawDescriptionHelpFormatter
20
21
22# define and zero variables used in the program here.
23totOxygens = 0
24totHydrogens = 0
25Hmat = []
26
27oxygenList = []
28
29
30
31def usage():
32 print(__doc__)
33
34class oxygenAtom:
35 def __init__(self):
36 self.pos_ = [0.0, 0.0, 0.0]
37 self.donors_ = []
38 self.acceptors_ = []
39 self.isSurface_ = 'false'
40 self.neighbors_ = []
41
42 def setPos(self, pos):
43 self.pos_ = pos
44
45 def getPos(self):
46 return self.pos_
47
48class proton:
49 def __init__(self):
50 self.pos_ = [0.0, 0.0, 0.0]
51
52 def setPos(self, pos):
53 self.pos_ = pos
54
55 def getPos(self):
56 return self.pos_
57
58
59def roundMe(x):
60 if (x >= 0.0):
61 return math.floor(x + 0.5)
62 else:
63 return math.ceil(x - 0.5)
64
65def wrapMe(disp, boxl):
66 scaled = 0.0
67 scaled = disp / boxl
68 scaled = scaled - roundMe(scaled)
69 disp = scaled * boxl
70 return disp
71
72def normalize(L1):
73 L2 = []
74 myLength = math.sqrt(dot(L1, L1))
75 for i in range(len(L1)):
76 L2.append(L1[i] / myLength)
77 return L2
78
79def dot(L1, L2):
80 myDot = 0.0
81 for i in range(len(L1)):
82 myDot = myDot + L1[i]*L2[i]
83 return myDot
84
85def addVectors(v1, v2):
86 newV = []
87 if (len(v1) != len(v2)):
88 print("Can't add vectors of different sizes!")
89 else:
90 for i in range(0, len(v1)):
91 newV.append(v1[i] + v2[i])
92 return newV
93
94def subtractVectors(v1, v2):
95 newV = []
96 if (len(v1) != len(v2)):
97 print("Can't subtract vectors of different sizes!")
98 else:
99 for i in range(0, len(v1)):
100 newV.append(v1[i] - v2[i])
101 return newV
102
103
104
105
106
107
108
109def readInputFile(inputFileName):
110 inputFile = open(inputFileName, "r")
111
112 line = inputFile.readline()
113 totOxygens = int(line.split()[0])
114 print("total number of oxygens to read in = ", totOxygens)
115
116 line = inputFile.readline()
117 Hxx = float(line.split()[1].strip(';'))
118 Hxy = float(line.split()[2].strip(';'))
119 Hxz = float(line.split()[3].strip(';'))
120 Hyx = float(line.split()[4].strip(';'))
121 Hyy = float(line.split()[5].strip(';'))
122 Hyz = float(line.split()[6].strip(';'))
123 Hzx = float(line.split()[7].strip(';'))
124 Hzy = float(line.split()[8].strip(';'))
125 Hzz = float(line.split()[9].strip(';'))
126 Hmat.append([Hxx, Hxy, Hxz])
127 Hmat.append([Hyx, Hyy, Hyz])
128 Hmat.append([Hzx, Hzy, Hzz])
129 print(Hmat)
130
131 while True:
132 line = inputFile.readline()
133 if not line:
134 break
135
136 newOxygen = oxygenAtom()
137 posVec = []
138 posX = float(line.split()[1])
139 posY = float(line.split()[2])
140 posZ = float(line.split()[3])
141 posVec.append([posX, posY, posZ])
142 newOxygen.setPos(posVec)
143
144 oxygenList.append(newOxygen)
145
146 print("total number of oxygens read in = ", len(oxygenList))
147 return (oxygenList, Hmat)
148
149
150
151
152def findNeighbors(oxygenList):
153 print("finding neighbors now")
154 for i in range(0, len(oxygenList)):
155 oxygenA = oxygenList[i]
156 for j in range(i+1, len(oxygenList)):
157 oxygenB = oxygenList[j]
158
159 if (oxygenA == oxygenB):
160 print("Querying same oxygen!")
161
162 # Want to wrap x & y displacements, but leave z as is
163 # since z exposes to a surface
164 xDisp = oxygenA.getPos()[0][0] - oxygenB.getPos()[0][0]
165 xDisp = wrapMe(xDisp, Hmat[0][0])
166 yDisp = oxygenA.getPos()[0][1] - oxygenB.getPos()[0][1]
167 yDisp = wrapMe(yDisp, Hmat[1][1])
168 zDisp = oxygenA.getPos()[0][2] - oxygenB.getPos()[0][2]
169 #zDisp = wrapMe(zDisp,Hmat[2][2])
170
171 dist = np.sqrt( pow(xDisp, 2) + pow(yDisp, 2) + pow(zDisp, 2) )
172
173 if (dist < 3.0):
174 oxygenA.neighbors_.append(j)
175 oxygenB.neighbors_.append(i)
176
177 nSurfaceAtoms = 0
178 for i in range(0, len(oxygenList)):
179 oxygenA = oxygenList[i]
180 if ( len(oxygenA.neighbors_) < 4):
181 nSurfaceAtoms = nSurfaceAtoms + 1
182 oxygenA.isSurface_ = 'true'
183 oxygenA.neighbors_.append(-1)
184 print("nSurfaceAtoms =", nSurfaceAtoms)
185
186 #sanity check to make sure no-one is a self-neighbor
187 for i in range(0, len(oxygenList)):
188 oxygenA = oxygenList[i]
189 if i in oxygenA.neighbors_ :
190 print(i, oxygenA.neighbors_)
191
192
193def randomizeProtons(oxygenList):
194
195 print("Randomizing protons, begin bucket brigades")
196 # select random oxygen to start with and perform protonBucketBrigades.
197 for i in range(0, 1000):
198 startingOxygenIndex = random.randint(0, len(oxygenList)-1)
199 protonBucketBrigade(oxygenList, startingOxygenIndex)
200 if (i%100 == 0):
201 print("attempted", i, "brigades")
202
203
204 # check to see if any molecules have to many or still need protons
205 bulkUnderCoord = []
206 surfUnderCoord = []
207 for i in range(0, len(oxygenList)):
208 oxygen = oxygenList[i]
209
210 if (len(oxygen.donors_) > 2):
211 pass
212 #print "oxygen no.", i, "has", len(oxygen.donors_),"donors"
213 elif (len(oxygen.acceptors_) > 2):
214 pass
215 #print "oxygen no.", i, "has", len(oxygen.acceptors_),"acceptors"
216
217 if (oxygen.isSurface_ == 'true'):
218 numProtons = len(oxygen.donors_) + len(oxygen.acceptors_)
219 if (numProtons < 3):
220 surfUnderCoord.append(i)
221 elif (oxygen.isSurface_ == 'false'):
222 numProtons = len(oxygen.donors_) + len(oxygen.acceptors_)
223 if (numProtons < 4):
224 bulkUnderCoord.append(i)
225
226 print("There are", len(bulkUnderCoord), "bulk oxygens needing protons")
227 print("There are", len(surfUnderCoord), "surf. oxygens needing protons")
228
229 print("Attempting more protonBucketBridages starting with the deficient oxygens")
230 print("starting with bulk undercoordinated oxygens")
231 for j in range(0, 10):
232 for i in range(0, len(bulkUnderCoord)):
233 protonBucketBrigade(oxygenList, bulkUnderCoord[i])
234 if (j%10 == 0):
235 print("attempted", j, "brigades")
236
237 print("moving onto surface undercoordinated oxygens")
238 for j in range(0, 10):
239 for i in range(0, len(surfUnderCoord)):
240 protonBucketBrigade(oxygenList, surfUnderCoord[i])
241 if (j%10 == 0):
242 print("attempted", j, "brigades")
243
244
245 # re-check to see if any molecules have to many or still need protons
246 bulkUnderCoord = []
247 surfUnderCoord = []
248 for i in range(0, len(oxygenList)):
249 oxygen = oxygenList[i]
250
251 if (oxygen.isSurface_ == 'true'):
252 numProtons = len(oxygen.donors_) + len(oxygen.acceptors_)
253 if (numProtons < 3):
254 surfUnderCoord.append(i)
255 print('surfUC', i, oxygen.neighbors_, oxygen.donors_, oxygen.acceptors_)
256 elif (oxygen.isSurface_ == 'false'):
257 numProtons = len(oxygen.donors_) + len(oxygen.acceptors_)
258 if (numProtons < 4):
259 bulkUnderCoord.append(i)
260 print('bulkUC', i, oxygen.neighbors_, oxygen.donors_, oxygen.acceptors_)
261
262
263 print("There are", len(bulkUnderCoord), "bulk oxygens needing protons")
264 print("There are", len(surfUnderCoord), "surf. oxygens needing protons")
265
266
267 #Let's check to see the total number of donors_ and acceptors_ cites on the surfaces.
268 # If these numbers aren't equal, the lattice can end up frustrated (lacking or excess of protons)
269 numDonors = 0
270 numAcceptors = 0
271 for i in range(0, len(oxygenList)):
272 oxygen = oxygenList[i]
273 if (oxygen.isSurface_ == 'true'):
274 numDonors = numDonors + len(oxygen.donors_)
275 numAcceptors = numAcceptors + len(oxygen.acceptors_)
276 print("numDonors =", numDonors)
277 print("numAcceptors =", numAcceptors)
278
279 for i in range(0, len(oxygenList)):
280 oxygen = oxygenList[i]
281
282 if (oxygen.isSurface == 'true'):
283 numBonds = len(oxygen.donors_) + len(oxygen.acceptors_)
284 if (numBonds < 3):
285 print('\nSurfDeficient')
286 print('UCoxygen', i, oxygen.neighbors_, oxygen.donors_, oxygen.acceptors_)
287 for i in range(0, len(oxygen.neighbors_)):
288 nborOx = oxygenList[oxygen.neighbors_[i]]
289 print("neighbor", oxygen.neighbors_[i], nborOx.neighbors_, nborOx.donors_, nborOx.acceptors_)
290
291 elif (oxygen.isSurface_ == 'false'):
292 if (len(oxygen.donors_) < 2):
293 print('\nDonorDeficient')
294 print("UCoxygen", i, oxygen.neighbors_, oxygen.donors_, oxygen.acceptors_)
295 for i in range(0, len(oxygen.neighbors_)):
296 nborOx = oxygenList[oxygen.neighbors_[i]]
297 print("neighbor", oxygen.neighbors_[i], nborOx.neighbors_, nborOx.donors_, nborOx.acceptors_)
298 if (len(oxygen.acceptors_) < 2):
299 print('\nAcceptorDeficient')
300 print("UCoxygen", i, oxygen.neighbors_, oxygen.donors_, oxygen.acceptors_)
301 for i in range(0, len(oxygen.neighbors_)):
302 nborOx = oxygenList[oxygen.neighbors_[i]]
303 print("neighbor", oxygen.neighbors_[i], nborOx.neighbors_, nborOx.donors_, nborOx.acceptors_)
304
305
306
307def protonBucketBrigade(oxygenList, startOxygenIndex):
308 # First we need to loop through and create a dummy set of the oxygen. These will be tested and incremented for Hbonds
309 dummyList = []
310 for i in range(0, len(oxygenList)):
311 dummyOxygen = copy.deepcopy(oxygenList[i])
312 dummyList.append(dummyOxygen)
313
314 # assign starting oxygen to a dummy starting oxygen
315 dummyStart = dummyList[startOxygenIndex]
316 dummyStartIndex = startOxygenIndex
317
318
319 # badNeighbors are oxygens already donating to/accepting from current
320 badNbors = set(dummyStart.donors_).union(set(dummyStart.acceptors_))
321 # goodNbors are my nbors not in badnbors
322 goodNbors = set(dummyStart.neighbors_).difference(badNbors)
323
324
325 if (len(goodNbors) == 0):
326 return
327
328 # sample next oxygen from good neighbors
329 if (dummyStart.isSurface_ == 'true'):
330 numBonds = len(dummyStart.donors_) + len(dummyStart.acceptors_)
331 if (numBonds < 3):
332 if ( len(dummyStart.donors_) < 2):
333 dummyNextIndex = random.sample(goodNbors, 1)[0]
334 dummyNext = dummyList[dummyNextIndex]
335
336 if (dummyNext.isSurface_ == 'true'):
337 nOnumBonds = len(dummyNext.donors_) + len(dummyNext.acceptors_)
338 if (nOnumBonds < 3):
339 if (len(dummyNext.acceptors_) < 2):
340 dummyStart.donors_.append(dummyNextIndex)
341 dummyNext.acceptors_.append(dummyStartIndex)
342 elif (len(dummyNext.acceptors_) == 2):
343 return #nextOxygen has saturated acceptors
344 elif (nOnumBonds == 3):
345 return #nextOxygen has total saturation of bonds
346 elif (nOnumBonds > 3):
347 print("Code is broken, nextOxygen is surface and has > 3 bonds.")
348 print(dummyNextIndex, dummyNext.neighbors_. dummyNext.donors_, dummyNext.acceptors_)
349
350 elif (dummyNext.isSurface_ == 'false'):
351 if (len(dummyNext.acceptors_) < 2):
352 dummyStart.donors_.append(dummyNextIndex)
353 dummyNext.acceptors_.append(dummyStartIndex)
354 elif (len(dummyNext.acceptors_) == 2):
355 return
356 elif (len(dummyNext.acceptors_) > 2):
357 print("Code is broken, nextOxygen is bulk and has > 2 acceptors.")
358 print(dummyNextIndex, dummyNext.neighbors_. dummyNext.donors_, dummyNext.acceptors_)
359
360 # we are a surface oxygen, but we have 2 donors already.
361 elif ( len(dummyStart.donors_) == 2):
362 return
363
364 elif ( len(dummyStart.donors_) > 2):
365 print("Code is broken, startOxygen is surface and has > 2 donors")
366 print(dummyStartIndex, dummyStart.neighbors_. dummyStart.donors_, dummyStart.acceptors_)
367
368 # we are a surface oxygen, but have saturation of bonds.
369 elif (numBonds == 3):
370 return
371
372 elif (numBonds > 3):
373 print("Code is broken, startOxygen has numBonds > 3 ")
374 print(dummyStartIndex, dummyStart.neighbors_. dummyStart.donors_, dummyStart.acceptors_)
375
376 elif (dummyStart.isSurface_ == 'false'):
377 if (len(dummyStart.donors_) < 2): #startOxygen < 2 donors, therefore can make a bond
378 dummyNextIndex = random.sample(goodNbors, 1)[0]
379 dummyNext = dummyList[dummyNextIndex]
380
381 if (dummyNext.isSurface_ == 'true'):
382 nOnumBonds = len(dummyNext.donors_) + len(dummyNext.acceptors_)
383 if (nOnumBonds < 3):
384 if ( len(dummyNext.acceptors_) < 2): #requirement met, make a bond
385 dummyStart.donors_.append(dummyNextIndex)
386 dummyNext.acceptors_.append(dummyStartIndex)
387 elif ( len(dummyNext.acceptors_) == 2):
388 return # nextOxygen has 2 acceptors, no room for one more.
389 elif (nOnumBonds == 3):
390 return #nextOxygen has 3 bonds, no room for one more.
391 elif (nOnumBonds > 3):
392 print("Code is broken, nextOxygen is surface and has > 3 bonds.")
393 print(dummyNextIndex, dummyNext.neighbors_. dummyNext.donors_, dummyNext.acceptors_)
394
395 elif (dummyNext.isSurface_ == 'false'):
396 if (len(dummyNext.acceptors_) < 2): #requirement met, make a bond
397 dummyStart.donors_.append(dummyNextIndex)
398 dummyNext.acceptors_.append(dummyStartIndex)
399 elif (len(dummyNext.acceptors_) == 2):
400 return # nextOxygen has 2 acceptors, no room for one more.
401 elif (len(dummyNext.acceptors_) > 2):
402 print("Code is broken, nextOxygen is bulk and has > 2 acceptors.")
403 print(dummyNextIndex, dummyNext.neighbors_. dummyNext.donors_, dummyNext.acceptors_)
404
405
406 #startOxygen is not surface and has 2 donors, no more room for additional donors.
407 elif (len(startDummy.donors_) == 2):
408 return
409
410 elif (len(dummyStart.donors_) > 2):
411 print("Code broke, startOxygen is not surface and has > 2 donors.")
412 print(dummyStartIndex, dummyStart.neighbors_. dummyStart.donors_, dummyStart.acceptors_)
413
414
415 # update nextOxygen to currentOxygen, begin looping.
416 dummyPrev = dummyStart
417 dummyPrevIndex = dummyStartIndex
418 dummyCurrent = dummyNext
419 dummyCurrentIndex = dummyNextIndex
420
421
422 # Now we will loop until we reach back to the starting Oxygen, guaranteeing a zero-net-dipole chain!
423 while (dummyCurrentIndex != dummyStartIndex):
424
425 badNbors = set(dummyCurrent.donors_).union(set(dummyCurrent.acceptors_))
426 goodNbors = set(dummyCurrent.neighbors_).difference(badNbors)
427
428
429 # if currentOxygen is surface, must test previousOxygen for surface or bulk.
430 if (dummyCurrent.isSurface_ == 'true'):
431 #can't test the raw number of donors_ could have a (DAA) set of bonds and donors < 2 would give false positive.
432 numBonds = len(dummyCurrent.donors_) + len(dummyCurrent.acceptors_)
433 if (numBonds < 3): # with < 3 bonds, currentOxygen has room to donate to another oxygen.
434 dummyNextIndex = random.sample(goodNbors, 1)[0]
435 dummyNext = dummyList[dummyNextIndex]
436
437 # if previous and current both surface, nextOxygen must be a bulk.
438 if (dummyPrev.isSurface_ == 'true'):
439 if (dummyNext.isSurface_ == 'true'):
440 return
441 elif (dummyNext.isSurface_ == 'false'):
442 if (len(dummyNext.acceptors_) < 2):
443 dummyCurrent.donors_.append(dummyNextIndex)
444 dummyNext.acceptors_.append(dummyCurrentIndex)
445 elif (len(dummyNext.acceptors_) == 2):
446 return # nextOxygen has 2 acceptors, no room for one more.
447 elif (len(dummyNext.acceptors_) > 2):
448 print("Code is broken (while), nextOxygen is bulk and has > 2 acceptors.")
449
450
451 # if previousOxygen is not surface, then nextOxygen can be a surface or bulk oxygen
452 elif (dummyPrev.isSurface_ == 'false'):
453
454 if (dummyNext.isSurface_ == 'true'):
455 nOnumBonds = len(dummyNext.donors_) + len(dummyNext.acceptors_)
456 if (nOnumBonds < 3):
457 if (len(dummyNext.acceptors_) < 2):
458 dummyCurrent.donors_.append(dummyNextIndex)
459 dummyNext.acceptors_.append(dummyCurrentIndex)
460 elif (len(dummyNext.acceptors_) == 2):
461 return # nextOxygen has 2 acceptors, no room for one more.
462 elif (len(dummyNext.acceptors_) > 2):
463 print("Code broken (while), nextOxygen has > 2 acceptors")
464 elif (nOnumBonds == 3):
465 return #nextOxygen has 3 bonds, no room for one more.
466 elif (nOnumBonds > 3):
467 print("Code is broken (while), nextOxygen is surface and has > 3 bonds.")
468
469 elif (dummyNext.isSurface_ == 'false'):
470 if (len(dummyNext.acceptors_) < 2):
471 dummyCurrent.donors_.append(dummyNextIndex)
472 dummyNext.acceptors_.append(dummyCurrentIndex)
473 elif (len(dummyNext.acceptors_) == 2):
474 return # nextOxygen has 2 acceptors, no room for one more.
475 elif (len(dummyNext.acceptors_) > 2):
476 print("Code is broken (while), nextOxygen is bulk and has > 2 acceptors.")
477
478 elif (numBonds == 3):
479 return #currentOxygen has a saturation of bonds.
480 elif (numBonds > 3):
481 print("Code is broken (while), currentOxygen is surface and has > 3 bonds.")
482
483 #if currentOxygen is not surface, nextOxygen can be bulk or surface.
484 elif (dummyCurrent.isSurface_ == 'false'):
485 dummyNextIndex = random.sample(goodNbors, 1)[0]
486 dummyNext = dummyList[dummyNextIndex]
487
488
489 if (dummyNext.isSurface_ == 'true'):
490 nOnumBonds = len(dummyNext.donors_) + len(dummyNext.acceptors_)
491 if (nOnumBonds < 3):
492 if (len(dummyNext.acceptors_) < 2):
493 dummyCurrent.donors_.append(dummyNextIndex)
494 dummyNext.acceptors_.append(dummyCurrentIndex)
495 elif (len(dummyNext.acceptors_) == 2):
496 return # nextOxygen has 2 acceptors, no room for one more.
497 elif (len(dummyNext.acceptors_) > 2):
498 print("Code broken (while), nextOxygen has > 2 acceptors")
499 elif (nOnumBonds == 3):
500 return #nextOxygen has 3 bonds, no room for one more.
501 elif (nOnumBonds > 3):
502 print("Code is broken (while), nextOxygen is surface and has > 3 bonds.")
503
504 elif (dummyNext.isSurface_ == 'false'):
505 if (len(dummyNext.acceptors_) < 2):
506 dummyCurrent.donors_.append(dummyNextIndex)
507 dummyNext.acceptors_.append(dummyCurrentIndex)
508 elif (len(dummyNext.acceptors_) == 2):
509 return # nextOxygen has 2 acceptors, no room for one more.
510 elif (len(dummyNext.acceptors_) > 2):
511 print("Code is broken (while), nextOxygen is bulk and has > 2 acceptors.")
512
513
514 # update nextOxygen to currentOxygen, begin looping.
515 dummyPrev = dummyCurrent
516 dummyPrevIndex = dummyCurrentIndex
517 dummyCurrent = dummyNext
518 dummyCurrentIndex = dummyNextIndex
519
520 #post while loop: test for length of protonLoop, will be set to zero upon break above.
521 if (len(dummyList) > 0):
522 assignBonds(oxygenList, dummyList)
523
524
525
526
527def assignBonds(oxygenList, dummyList):
528 # If we are here, we have successfully made a chain of hydrogen-bonds.
529 # We will save this by overwriting the oxygenList with the dummyList
530 for i in range(0, len(dummyList)):
531 oxygenList[i] = copy.deepcopy(dummyList[i])
532
533
534
535
536
537def monteCarloMethod(oxygenList):
538
539 print("Assigning all oxygens a random set of bonds now")
540 #First we will assign all of the oxygens a random donor/acceptor set. We can end up with > 2 donors or acceptors here.
541 numProtons = 0
542 for i in range(0, len(oxygenList)):
543 oxygenA = oxygenList[i]
544
545 for j in range(0, len(oxygenA.neighbors_)):
546 oxygenBIndex = oxygenA.neighbors_[j]
547
548 if (oxygenBIndex not in oxygenA.donors_ and oxygenBIndex not in oxygenA.acceptors_):
549 rand = random.uniform(0, 1)
550 if (rand < 0.5):
551 oxygenA.donors_.append(oxygenA.neighbors_[j])
552 if (oxygenBIndex != -1):
553 oxygenList[oxygenA.neighbors_[j]].acceptors_.append(i)
554 if (rand >= 0.5):
555 oxygenA.acceptors_.append(oxygenA.neighbors_[j])
556 if (oxygenBIndex != -1):
557 oxygenList[oxygenA.neighbors_[j]].donors_.append(i)
558 numProtons = numProtons + len(oxygenA.donors_)
559 print("numProtons =", numProtons)
560
561 numUnOr = 0
562 for i in range(0, len(oxygenList)):
563 oxygen = oxygenList[i]
564 if (len(oxygen.donors_) < 2):
565 numUnOr = numUnOr + 1
566 print("Pre MC:", numUnOr, "oxygens with < 2 donors_")
567
568
569
570 print("Calculating net dipole of the system")
571 netDipole = [0.0, 0.0, 0.0]
572 #Next we need to calculate the dipole moment of the system
573 for i in range(0, len(oxygenList)):
574 oxygenA = oxygenList[i]
575
576 #donor configuration has dipole moment (+1) * rUnit vec between two oxygens
577 for j in range(0, len(oxygenA.donors_)):
578 oxygenB = oxygenList[oxygenA.donors_[j]]
579
580 rVec = subtractVectors(oxygenA.getPos()[0], oxygenB.getPos()[0])
581 rUnit = normalize(rVec)
582
583 netDipole = addVectors(netDipole, rUnit)
584
585 #acceptor configuration has dipole moment (-1) * rUnit vec between two oxygens
586 for j in range(0, len(oxygenA.acceptors_)):
587 oxygenB = oxygenList[oxygenA.acceptors_[j]]
588
589 rVec = subtractVectors(oxygenA.getPos()[0], oxygenB.getPos()[0])
590 rUnit = normalize(rVec)
591
592 netDipole = subtractVectors(netDipole, rUnit)
593
594 #Now we divide by 2.0 since we visited each pair twice
595 for i in range(0, len(netDipole)):
596 netDipole[i] = netDipole[i] / 2.0
597 print("net dipole pre MC moves =", netDipole)
598
599 totalCycles = 100000
600 nBarCycles = totalCycles / 5.0
601 for nCycles in range(0, totalCycles):
602 #update kT(nCycles)
603 kT = 10.0*np.exp(-nCycles/nBarCycles)
604
605
606 #Next we will perform MonteCarlo moves where we will flip the donors / acceptors
607 for i in range(0, len(oxygenList)):
608 #oxygenAIndex = random.randint(0,len(oxygenList)-1)
609 oxygenAIndex = i
610 oxygenA = oxygenList[oxygenAIndex]
611
612 oxygenBIndex = random.sample(oxygenA.neighbors_, 1)[0]
613 if (oxygenBIndex != -1):
614 oxygenB = oxygenList[oxygenBIndex]
615
616
617 #compute pre-flip dipole moment of Hbond
618 oldDipole = []
619
620 rVec = subtractVectors(oxygenA.getPos()[0], oxygenB.getPos()[0])
621 rUnit = normalize(rVec)
622
623 if (oxygenBIndex in oxygenA.donors_):
624 for i in range(0, len(netDipole)):
625 oldDipole.append(1.0 * rUnit[i])
626
627 elif (oxygenBIndex in oxygenA.acceptors_):
628 for i in range(0, len(netDipole)):
629 oldDipole.append(-1.0 * rUnit[i])
630
631
632 EpreFlip = calcPairEnergy(oxygenA, oxygenB, netDipole, kT)
633
634 swapDonorAcceptor(oxygenA, oxygenAIndex, oxygenB, oxygenBIndex)
635
636 #compute post-flip dipole moment of Hbond
637 newDipole = []
638
639 rVec = subtractVectors(oxygenA.getPos()[0], oxygenB.getPos()[0])
640 rUnit = normalize(rVec)
641
642 if (oxygenBIndex in oxygenA.donors_):
643 for i in range(0, len(netDipole)):
644 newDipole.append(1.0 * rUnit[i])
645
646 elif (oxygenBIndex in oxygenA.acceptors_):
647 for i in range(0, len(netDipole)):
648 newDipole.append(-1.0 * rUnit[i])
649
650
651 #compute new net dipole moment of system
652 newNetDipole = []
653 for i in range(0, len(netDipole)):
654 newNetDipole.append(netDipole[i] - oldDipole[i] + newDipole[i])
655
656 EpostFlip = calcPairEnergy(oxygenA, oxygenB, newNetDipole, kT)
657
658
659 #accept or reject move
660 deltaE = EpostFlip - EpreFlip
661 randNum = random.uniform(0, 1)
662
663 if (randNum <= np.exp(-deltaE / kT) ):
664 netDipole = newNetDipole
665 continue
666 else:
667 swapDonorAcceptor(oxygenA, oxygenAIndex, oxygenB, oxygenBIndex)
668
669 if (nCycles%100 == 0):
670 print(nCycles, kT, netDipole)
671
672
673 #after MC cycles, let's check to see how many oxygens have the right number of donors/acceptors
674 numUnOr = 0
675 for i in range(0, len(oxygenList)):
676 oxygen = oxygenList[i]
677 if (len(oxygen.donors_) < 2 and len(oxygen.acceptors_) < 2):
678 numUnOr = numUnOr + 1
679 print("Post MC:", numUnOr, "oxygens with < 2 donors and < 2 acceptors")
680
681
682
683def calcPairEnergy(oxygenA, oxygenB, netDipole, kT):
684 #f = weighting for having 2 donors & 2 acceptors
685 #g = weighting for zero-net dipole to the system
686 f = 2.0/kT
687 g = 1.0/kT
688
689 oxAda = pow(len(oxygenA.donors_) - len(oxygenA.acceptors_), 2.0)
690 oxBda = pow(len(oxygenB.donors_) - len(oxygenB.acceptors_), 2.0)
691
692 E = f*( oxAda + oxBda ) + g*(np.sqrt(dot(netDipole, netDipole)))
693
694 return E
695
696
697def swapDonorAcceptor(oxygenA, oxygenAIndex, oxygenB, oxygenBIndex):
698 if (oxygenBIndex in oxygenA.donors_):
699
700 listAIndex = oxygenA.donors_.index(oxygenBIndex)
701 listBIndex = oxygenB.acceptors_.index(oxygenAIndex)
702
703 del oxygenA.donors_[listAIndex]
704 del oxygenB.acceptors_[listBIndex]
705
706 oxygenA.acceptors_.append(oxygenBIndex)
707 oxygenB.donors_.append(oxygenAIndex)
708
709 elif (oxygenBIndex in oxygenA.acceptors_):
710
711 listAIndex = oxygenA.acceptors_.index(oxygenBIndex)
712 listBIndex = oxygenB.donors_.index(oxygenAIndex)
713
714 del oxygenA.acceptors_[listAIndex]
715 del oxygenB.donors_[listBIndex]
716
717 oxygenA.donors_.append(oxygenBIndex)
718 oxygenB.acceptors_.append(oxygenAIndex)
719
720
721
722
723def addProtonsToDonors(oxygenList):
724 print("Adding protons to the lattice")
725
726 numProtons = 0
727 protonList = []
728 for i in range(0, len(oxygenList)):
729 oxygenA = oxygenList[i]
730 numProtons = numProtons + len(oxygenA.donors_)
731 #Should be able to just scan all the donors_, if someone is their second donor, then another oxygen must have it.
732 for j in range(0, len(oxygenA.donors_)):
733 oxygenBIndex = oxygenA.donors_[j]
734
735 if (oxygenBIndex != -1):
736 oxygenB = oxygenList[oxygenA.donors_[j]]
737
738 # place protons 1 unit vector away from the donating oxygen
739 # towards the accepting oxygen
740 xDisp = oxygenA.getPos()[0][0] - oxygenB.getPos()[0][0]
741 yDisp = oxygenA.getPos()[0][1] - oxygenB.getPos()[0][1]
742 zDisp = oxygenA.getPos()[0][2] - oxygenB.getPos()[0][2]
743
744 rVec = [xDisp, yDisp, zDisp]
745 rUnit = normalize(rVec)
746
747 xPos = oxygenA.getPos()[0][0] + rUnit[0]
748 yPos = oxygenA.getPos()[0][1] + rUnit[1]
749 zPos = oxygenA.getPos()[0][2] + rUnit[2]
750 protonVec = [xPos, yPos, zPos]
751
752 newProton = proton()
753 newProton.setPos(protonVec)
754
755 protonList.append(newProton)
756
757 elif (oxygenBIndex == -1):
758 # for the (-1) neigbors, we have to take the negatice of the average of the other 3 HBond vectors
759 aveVec = [0.0, 0.0, 0.0]
760 for k in range(0, len(oxygenA.neighbors_)):
761 oxNeighIndex = oxygenA.neighbors_[k]
762
763 if (oxNeighIndex != -1):
764 oxygenB = oxygenList[oxygenA.donors_[j]]
765
766 xDisp = oxygenA.getPos()[0][0] - oxygenB.getPos()[0][0]
767 yDisp = oxygenA.getPos()[0][1] - oxygenB.getPos()[0][1]
768 zDisp = oxygenA.getPos()[0][2] - oxygenB.getPos()[0][2]
769
770 rVec = [xDisp, yDisp, zDisp]
771 aveVec = addVectors(aveVec, rVec)
772
773 #normalize aveVec, and multiply by (-1)
774 aveVec = normalize(aveVec)
775 for k in range(0, len(aveVec)):
776 aveVec[k] = -1.0 * aveVec[k]
777
778 xPos = oxygenA.getPos()[0][0] + aveVec[0]
779 yPos = oxygenA.getPos()[0][1] + aveVec[1]
780 zPos = oxygenA.getPos()[0][2] + aveVec[2]
781 protonVec = [xPos, yPos, zPos]
782
783 newProton = proton()
784 newProton.setPos(protonVec)
785
786 protonList.append(newProton)
787
788 print("numProtons = ", numProtons)
789 print("len(oxygenList) =", len(oxygenList))
790 print("len(protonList) =", len(protonList))
791 return protonList
792
793
794
795def computeQuats(oxygenList):
796
797 #for i in range(len(indices)):
798 # DonatingTo.append(list())
799 # AcceptingFrom.append(list())
800 # OldDonor.append(list())
801 # print "Computing Quaternions"
802
803 # Set up initial rotation matrix
804 ux = [0.0, 0.0, 0.0]
805 uy = [0.0, 0.0, 0.0]
806 uz = [0.0, 0.0, 0.0]
807 RotMat = [ux, uy, uz]
808 totalDipole = [0.0, 0.0, 0.0]
809 #for i in range(len(indices)):
810 # print "doing quats for molecule %d" % i
811 # print "Dlen = %d " % len(DonatingTo[i])
812 # print DonatingTo[i]
813 for i in range(0, len(oxygenList)):
814 oxygenA = oxygenList[i]
815 #print oxygenA.donors_
816 #print oxygenA.acceptors_
817
818 #myPos = positions[i]
819
820 #acceptor1 = DonatingTo[i][0]
821 #acceptor2 = DonatingTo[i][1]
822
823 #acceptor1 = oxygenA.donors_[0]
824 #acceptor2 = oxygenA.donors_[1]
825
826
827 '''
828 if (acceptor1 == -1 and acceptor2 == -1):
829 donor1 = AcceptingFrom[i][0]
830 donor2 = AcceptingFrom[i][1]
831 npos = positions[donor1]
832 apos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
833 apos1 = wrapVector(apos1)
834 npos = positions[donor2]
835 apos2 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
836 apos2 = wrapVector(apos2)
837 tempX = [0.0, 0.0, 0.0]
838 tempY = [0.0, 0.0, 0.0]
839 tempZ = [0.0, 0.0, 0.0]
840 tempZ[0] = 0.5*(apos1[0] + apos2[0])
841 tempZ[1] = 0.5*(apos1[1] + apos2[1])
842 tempZ[2] = 0.5*(apos1[2] + apos2[2])
843 tempX[0] = (apos2[0] - apos1[0])
844 tempX[1] = (apos2[1] - apos1[1])
845 tempX[2] = (apos2[2] - apos1[2])
846 tempZ = normalize(tempZ)
847 tempX = normalize(tempX)
848 tempY = cross(tempX, tempZ)
849 dpos1[0] = tempZ[0] - 0.5*tempY[0]
850 dpos1[1] = tempZ[1] - 0.5*tempY[1]
851 dpos1[2] = tempZ[2] - 0.5*tempY[2]
852 dpos2[0] = tempZ[0] + 0.5*tempY[0]
853 dpos2[1] = tempZ[1] + 0.5*tempY[1]
854 dpos2[2] = tempZ[2] + 0.5*tempY[2]
855
856 else:
857 if (acceptor1 == -1):
858 tempVec = [0.0, 0.0, 0.0]
859 for j in range(3):
860 thisNeighbor = neighbors[i][j]
861 npos = positions[thisNeighbor]
862 npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
863 npos1 = wrapVector(npos1)
864 tempVec[0] = tempVec[0] + npos1[0]
865 tempVec[1] = tempVec[1] + npos1[1]
866 tempVec[2] = tempVec[2] + npos1[2]
867 dpos1 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
868 dpos1 = normalize(dpos1)
869 else:
870 a1pos = positions[acceptor1]
871 dpos1 = [a1pos[0] - myPos[0], a1pos[1] - myPos[1], a1pos[2] - myPos[2]]
872 dpos1 = wrapVector(dpos1)
873 dpos1 = normalize(dpos1)
874
875
876 if (acceptor2 == -1):
877 tempVec = [0.0, 0.0, 0.0]
878 for j in range(3):
879 thisNeighbor = neighbors[i][j]
880 npos = positions[thisNeighbor]
881 npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
882 npos1 = wrapVector(npos1)
883 tempVec[0] = tempVec[0] + npos1[0]
884 tempVec[1] = tempVec[1] + npos1[1]
885 tempVec[2] = tempVec[2] + npos1[2]
886 dpos2 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
887 dpos2 = normalize(dpos2)
888 else:
889 a2pos = positions[acceptor2]
890 dpos2 = [a2pos[0] - myPos[0], a2pos[1] - myPos[1], a2pos[2] - myPos[2]]
891 dpos2 = wrapVector(dpos2)
892 dpos2 = normalize(dpos2)
893
894 for j in range(3):
895 uz[j] = (dpos1[j] + dpos2[j])/2.0
896 uz = normalize(uz)
897 for j in range(3):
898 uy[j] = dpos2[j] - dpos1[j]
899 uy = normalize(uy)
900 ux = cross(uy, uz)
901 ux = normalize(ux)
902
903 q = [0.0, 0.0, 0.0, 0.0]
904
905 # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
906
907 RotMat[0] = ux
908 RotMat[1] = uy
909 RotMat[2] = uz
910
911 t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
912
913 if( t > 1e-6 ):
914 s = 0.5 / math.sqrt( t )
915 q[0] = 0.25 / s
916 q[1] = (RotMat[1][2] - RotMat[2][1]) * s
917 q[2] = (RotMat[2][0] - RotMat[0][2]) * s
918 q[3] = (RotMat[0][1] - RotMat[1][0]) * s
919 else:
920 ad1 = RotMat[0][0]
921 ad2 = RotMat[1][1]
922 ad3 = RotMat[2][2]
923
924 if( ad1 >= ad2 and ad1 >= ad3 ):
925 s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
926 q[0] = (RotMat[1][2] - RotMat[2][1]) * s
927 q[1] = 0.25 / s
928 q[2] = (RotMat[0][1] + RotMat[1][0]) * s
929 q[3] = (RotMat[0][2] + RotMat[2][0]) * s
930 elif ( ad2 >= ad1 and ad2 >= ad3 ):
931 s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
932 q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
933 q[1] = (RotMat[0][1] + RotMat[1][0]) * s
934 q[2] = 0.25 / s
935 q[3] = (RotMat[1][2] + RotMat[2][1]) * s
936 else:
937 s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
938 q[0] = (RotMat[0][1] - RotMat[1][0]) * s
939 q[1] = (RotMat[0][2] + RotMat[2][0]) * s
940 q[2] = (RotMat[1][2] + RotMat[2][1]) * s
941 q[3] = 0.25 / s
942
943 quaternions[i] = q
944 totalDipole = [totalDipole[0] + uz[0], totalDipole[1] + uz[1], totalDipole[2] + uz[2]]
945 totalDipole = [-totalDipole[0], -totalDipole[1], -totalDipole[2]]
946 print totalDipole
947 Dipole = math.sqrt(dot(totalDipole, totalDipole))
948 print 'Total Dipole Moment = %d' % Dipole
949 if (Dipole > 40 or Dipole < -40):
950 print "Bad Dipole, starting over"
951 for i in range(len(indices)):
952 del OldDonor[:]
953 del AcceptingFrom[:]
954 del DonatingTo[:]
955 # del badChoices[:]
956 randomizeProtons()
957 computeQuats()
958 else:
959 print "All Done!"
960
961 '''
962
963
964
965def writeXYZfile(oxygenList, protonList, Hmat):
966 outFile = open("test.xyz", "w")
967
968 numObjects = len(oxygenList) + len(protonList)
969 outFile.write(str(numObjects) + "\n")
970 outFile.write("\t%f%s\t%f\t%f\t%f%s\t%f\t%f\t%f%s\t%f\t%f\t%f\n" %
971 (0.00000, ";",
972 Hmat[0][0], Hmat[0][1], Hmat[0][2], ";",
973 Hmat[1][0], Hmat[1][1], Hmat[1][2], ";",
974 Hmat[2][0], Hmat[2][1], Hmat[2][2]))
975 for i in range(0, len(oxygenList)):
976 oxygenA = oxygenList[i]
977 outFile.write("%s\t%f\t%f\t%f\n" % ('O', oxygenA.getPos()[0][0], oxygenA.getPos()[0][1], oxygenA.getPos()[0][2]))
978
979 for i in range(0, len(protonList)):
980 protonA = protonList[i]
981 outFile.write("%s\t%f\t%f\t%f\n" % ('H', protonA.getPos()[0], protonA.getPos()[1], protonA.getPos()[2]))
982
983
984# This function writes out the oxygens with their neighbors, donors, and acceptors
985def writeOxygenfile(oxygenList):
986 outFile = open("oxygen.dat", "w")
987
988 for i in range(0, len(oxygenList)):
989 oxygenA = oxygenList[i]
990 outFile.write(str(i) + "\t" + str(oxygenA.neighbors_) + "\t" + str(oxygenA.donors_) + "\t" + str(oxygenA.acceptors_) + "\n")
991
992
993def main(argv):
994 parser = argparse.ArgumentParser(
995 description='OpenMD module that generates proton disordered ice-Ih lattices given an (.xyz) coordinate file for the Oxygen lattice.',
996 formatter_class=RawDescriptionHelpFormatter,
997 epilog="Example: protonSampler -i oxygenLattice.xyz -o protonDisorderedIce.omd")
998 parser.add_argument("-i", "--inputFile=", action="store",
999 dest="inputFileName", help="use specified input file")
1000 parser.add_argument("-o", "--outputfile=", action="store", dest="outputFileName", help="use specified output (.omd) file")
1001
1002 if len(sys.argv) == 1:
1003 parser.print_help()
1004 sys.exit(2)
1005 args = parser.parse_args()
1006
1007 if (not args.inputFileName):
1008 parser.error("No inputFile was specified")
1009
1010 if (not args.outputFileName):
1011 parser.print_help()
1012 parser.error("No outputFile was specified")
1013
1014
1015
1016 (oxygenList, Hmat) = readInputFile(args.inputFileName)
1017 findNeighbors(oxygenList)
1018
1019 #randomizeProtons(oxygenList)
1020 monteCarloMethod(oxygenList)
1021
1022 protonList = addProtonsToDonors(oxygenList)
1023 writeXYZfile(oxygenList, protonList, Hmat)
1024 writeOxygenfile(oxygenList)
1025
1026 computeQuats(oxygenList)
1027
1028
1029 writeXYZfile(oxygenList, protonList, Hmat)
1030
1031
1032if __name__ == "__main__":
1033 #if len(sys.argv) == 1:
1034 #parser.print_help()
1035 #usage() # need to change to call argeparse stuffs
1036 #sys.exit()
1037 main(sys.argv[1:])