ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/utilities/waterRotator
(Generate patch)

Comparing trunk/src/applications/utilities/waterRotator (file contents):
Revision 1369 by cpuglis, Mon Nov 26 19:23:36 2007 UTC vs.
Revision 1370 by gezelter, Mon Oct 19 14:29:14 2009 UTC

# Line 18 | Line 18 | __author__ = "Dan Gezelter (gezelter@nd.edu)"
18   """
19  
20   __author__ = "Dan Gezelter (gezelter@nd.edu)"
21 < __version__ = "$Revision: 1.1 $"
22 < __date__ = "$Date: 2007-11-26 19:23:36 $"
21 > __version__ = "$Revision: 1.2 $"
22 > __date__ = "$Date: 2009-10-19 14:29:14 $"
23   __copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
24   __license__ = "OOPSE"
25  
# Line 28 | Line 28 | import random
28   import string
29   import math
30   import random
31 – from sets import *
32 – #from Numeric import *
31  
32   _haveMDFileName = 0
33   _haveOutputFileName = 0
# Line 54 | Line 52 | availableneighbs = []
52   DipoleVects = []
53   allAcceptors = []
54   availableneighbs = []
55 < surfaceSet = Set([-1])
55 > surfaceSet = set([-1])
56  
57  
58  
# Line 247 | Line 245 | def findNeighbors():
245  
246      surfaceCount = 0
247      for i in range(len(indices)):
248 <        if (len(neighbors[i]) == 3):
248 >        if (len(neighbors[i]) < 4):
249              neighbors[i].append(-1)
250              surfaceCount = surfaceCount + 1
251  
# Line 287 | Line 285 | def protonBucketBrigade(origPoint):
285          AcceptingFrom.append(list())
286          OldDonor.append(list())
287  
290 –    # print "randomizing Quaternions"
291 –    nIdeal = len(idealQuats)
292 –
288      donor = origPoint
289 <    # pick unreasonable start values (that don't match surface atom unreasonable values)
289 >    # pick unreasonable start values (that don't match surface atom
290 >    # unreasonable values)
291      acceptor = -10
296 –    #oldDonor = -10
297 –    # print 'origPoint = %d' % (origPoint)
292  
293      for i in range(50000):
294      #while (acceptor != origPoint):
295 <        myNeighbors = Set(neighbors[donor])
295 >        myNeighbors = set(neighbors[donor])
296              
297          # can't pick a proton choice from one of my current protons:
298 <        badChoices = Set(DonatingTo[donor]).union(Set(OldDonor[acceptor]))
298 >        badChoices = set(DonatingTo[donor]).union(set(OldDonor[acceptor]))
299          # can't send a proton back to guy who sent one to me (no give-backs):
300 <        #badChoices.add(Set(OldDonor[acceptor]))
300 >        #badChoices.add(set(OldDonor[acceptor]))
301          # can't send a proton to anyone who is already taking 2:
302          for j in myNeighbors.difference(surfaceSet):
303              if (len(AcceptingFrom[j]) == 2):
# Line 312 | Line 306 | def protonBucketBrigade(origPoint):
306          nDonors = len(DonatingTo[donor])
307  
308          if (nDonors <= 1):
309 <            acceptor = random.choice(list( myNeighbors.difference(badChoices) ) )
309 >
310 >            if (len(myNeighbors.difference(badChoices)) != 0):
311 >                acceptor = random.choice(list(myNeighbors.difference(badChoices)))
312 >            else:
313 >                acceptor = -1
314 >
315              DonatingTo[donor].append(acceptor)
316 +
317              if (acceptor != -1):
318                  AcceptingFrom[acceptor].append(donor)
319                  OldDonor[acceptor].append(donor)
# Line 322 | Line 322 | def protonBucketBrigade(origPoint):
322          else:
323              print "Whoah!  How'd we get here?"
324  
325 –        #OldDonor = donor
325          donor = acceptor
326          if (acceptor == -1):
327              # surface atoms all have a -1 neighbor, but that's OK.  A proton
# Line 353 | Line 352 | def computeQuats():
352          myPos = positions[i]
353  
354          acceptor1 = DonatingTo[i][0]
356 –        if (acceptor1 == -1):
357 –            tempVec = [0.0, 0.0, 0.0]
358 –            for j in range(3):
359 –                thisNeighbor = neighbors[i][j]
360 –                npos = positions[thisNeighbor]
361 –                npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
362 –                npos1 = wrapVector(npos1)
363 –                tempVec[0] = tempVec[0]  + npos1[0]
364 –                tempVec[1] = tempVec[1]  + npos1[1]
365 –                tempVec[2] = tempVec[2]  + npos1[2]                
366 –            dpos1 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
367 –            dpos1 = normalize(dpos1)
368 –        else:
369 –            a1pos = positions[acceptor1]
370 –            dpos1 = [a1pos[0] - myPos[0], a1pos[1] - myPos[1], a1pos[2] - myPos[2]]
371 –            dpos1 = wrapVector(dpos1)
372 –            dpos1 = normalize(dpos1)
373 –        
355          acceptor2 = DonatingTo[i][1]
356 <        if (acceptor2 == -1):
357 <            tempVec = [0.0, 0.0, 0.0]
358 <            for j in range(3):
359 <                thisNeighbor = neighbors[i][j]
360 <                npos = positions[thisNeighbor]
361 <                npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
362 <                npos1 = wrapVector(npos1)
363 <                tempVec[0] = tempVec[0]  + npos1[0]
364 <                tempVec[1] = tempVec[1]  + npos1[1]
365 <                tempVec[2] = tempVec[2]  + npos1[2]                
366 <            dpos2 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
367 <            dpos2 = normalize(dpos2)
368 <        else:
369 <            a2pos = positions[acceptor2]
370 <            dpos2 = [a2pos[0] - myPos[0], a2pos[1] - myPos[1], a2pos[2] - myPos[2]]
371 <            dpos2 = wrapVector(dpos2)
372 <            dpos2 = normalize(dpos2)
373 <
374 <
375 <
356 >        if (acceptor1 == -1 and acceptor2 == -1):
357 >            donor1 = AcceptingFrom[i][0]            
358 >            donor2 = AcceptingFrom[i][1]
359 >            npos = positions[donor1]
360 >            apos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
361 >            apos1 = wrapVector(apos1)
362 >            npos = positions[donor2]
363 >            apos2 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
364 >            apos2 = wrapVector(apos2)
365 >            tempX = [0.0, 0.0, 0.0]
366 >            tempY = [0.0, 0.0, 0.0]
367 >            tempZ = [0.0, 0.0, 0.0]
368 >            tempZ[0] = 0.5*(apos1[0] + apos2[0])
369 >            tempZ[1] = 0.5*(apos1[1] + apos2[1])
370 >            tempZ[2] = 0.5*(apos1[2] + apos2[2])
371 >            tempX[0] = (apos2[0] - apos1[0])
372 >            tempX[1] = (apos2[1] - apos1[1])
373 >            tempX[2] = (apos2[2] - apos1[2])
374 >            tempZ = normalize(tempZ)
375 >            tempX = normalize(tempX)
376 >            tempY = cross(tempX, tempZ)
377 >            dpos1[0] = tempZ[0] - 0.5*tempY[0]
378 >            dpos1[1] = tempZ[1] - 0.5*tempY[1]
379 >            dpos1[2] = tempZ[2] - 0.5*tempY[2]
380 >            dpos2[0] = tempZ[0] + 0.5*tempY[0]
381 >            dpos2[1] = tempZ[1] + 0.5*tempY[1]
382 >            dpos2[2] = tempZ[2] + 0.5*tempY[2]
383 >            
384 >        else:
385 >            if (acceptor1 == -1):
386 >                tempVec = [0.0, 0.0, 0.0]
387 >                for j in range(3):
388 >                    thisNeighbor = neighbors[i][j]
389 >                    npos = positions[thisNeighbor]
390 >                    npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
391 >                    npos1 = wrapVector(npos1)
392 >                    tempVec[0] = tempVec[0]  + npos1[0]
393 >                    tempVec[1] = tempVec[1]  + npos1[1]
394 >                    tempVec[2] = tempVec[2]  + npos1[2]                
395 >                dpos1 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
396 >                dpos1 = normalize(dpos1)
397 >            else:
398 >                a1pos = positions[acceptor1]
399 >                dpos1 = [a1pos[0] - myPos[0], a1pos[1] - myPos[1], a1pos[2] - myPos[2]]
400 >                dpos1 = wrapVector(dpos1)
401 >                dpos1 = normalize(dpos1)
402          
403  
404 +            if (acceptor2 == -1):
405 +                tempVec = [0.0, 0.0, 0.0]
406 +                for j in range(3):
407 +                    thisNeighbor = neighbors[i][j]
408 +                    npos = positions[thisNeighbor]
409 +                    npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
410 +                    npos1 = wrapVector(npos1)
411 +                    tempVec[0] = tempVec[0]  + npos1[0]
412 +                    tempVec[1] = tempVec[1]  + npos1[1]
413 +                    tempVec[2] = tempVec[2]  + npos1[2]                
414 +                dpos2 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
415 +                dpos2 = normalize(dpos2)
416 +            else:
417 +                a2pos = positions[acceptor2]
418 +                dpos2 = [a2pos[0] - myPos[0], a2pos[1] - myPos[1], a2pos[2] - myPos[2]]
419 +                dpos2 = wrapVector(dpos2)
420 +                dpos2 = normalize(dpos2)        
421 +
422          for j in range(3):
423              uz[j] = (dpos1[j] + dpos2[j])/2.0
424          uz = normalize(uz)
# Line 449 | Line 474 | def computeQuats():
474      print totalDipole
475      Dipole = math.sqrt(dot(totalDipole, totalDipole))
476      print 'Total Dipole Moment = %d' % Dipole
477 <    if (Dipole > 20 or Dipole < -20):
477 >    if (Dipole > 40 or Dipole < -40):
478          print "Bad Dipole, starting over"
479          for i in range(len(indices)):
480              del OldDonor[:]

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines