[ViewVC] Annotation of: OpenMD/branches/development/src/applications/utilities/waterRotator

#!/usr/bin/env python
"""Water Quaternion Sampler

Samples water orientations from a list of known good
orientations

Usage: waterRotator

Options:
  -h, --help              show this help
  -m, --meta-data=...     use specified meta-data (.md) file
  -o, --output-file=...   use specified output (.md) file


Example:
   waterRotator -m basal.md -o basal.new.md

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.1 $"
__date__ = "$Date: 2007-11-26 19:23:36 $"
__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
__license__ = "OOPSE"

import sys
import getopt
import string
import math
import random
from sets import *
#from Numeric import *

_haveMDFileName = 0
_haveOutputFileName = 0

metaData = []
frameData = []
positions = []
velocities = []
quaternions = []
angVels = []
indices = []
neighbors = []
DonatingTo = []
AcceptingFrom = []
OldDonor = []
Hmat = []
BoxInv = []
#Hmat = zeros([3,3],Float)
#BoxInv = zeros([3],Float)
H1vects = []
H2vects = []
DipoleVects = []
allAcceptors = []
availableneighbs = []
surfaceSet = Set([-1])


def usage():
    print __doc__


def readFile(mdFileName):
    mdFile = open(mdFileName, 'r')        
    # Find OOPSE version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line:
            OOPSEversion = line
            break
        line = mdFile.readline()
        
        # Rewind file and find start of MetaData block
        
    mdFile.seek(0)
    line = mdFile.readline()
    print "reading MetaData"
    while 1:
        if '<MetaData>' in line:
            while 2:
                metaData.append(line)
                line = mdFile.readline()
                if '</MetaData>' in line:
                    metaData.append(line)
                    break
            break
        line = mdFile.readline()

    mdFile.seek(0)
    print "reading Snapshot"
    line = mdFile.readline()
    while 1:
        if '<Snapshot>' in line:
            line = mdFile.readline()
            while 1:
                print "reading FrameData"
                if '<FrameData>' in line:
                    while 2:
                        frameData.append(line)
                        if 'Hmat:' in line:
                            L = line.split()
                            Hxx = float(L[2].strip(','))
                            Hxy = float(L[3].strip(','))
                            Hxz = float(L[4].strip(','))
                            Hyx = float(L[7].strip(','))
                            Hyy = float(L[8].strip(','))
                            Hyz = float(L[9].strip(','))
                            Hzx = float(L[12].strip(','))
                            Hzy = float(L[13].strip(','))
                            Hzz = float(L[14].strip(','))
                            Hmat.append([Hxx, Hxy, Hxz])
                            Hmat.append([Hyx, Hyy, Hyz])
                            Hmat.append([Hzx, Hzy, Hzz])
                            print Hmat
                            BoxInv.append(1.0/Hxx)
                            BoxInv.append(1.0/Hyy)
                            BoxInv.append(1.0/Hzz)
                            print BoxInv
                        line = mdFile.readline()
                        if '</FrameData>' in line:
                            frameData.append(line)
                            break
                    break

            line = mdFile.readline()
            while 1:
                if '<StuntDoubles>' in line:
                    line = mdFile.readline()
                    while 2:
                        L = line.split()
                        myIndex = int(L[0])
                        indices.append(myIndex)
                        pvqj = L[1]
                        x = float(L[2])
                        y = float(L[3])
                        z = float(L[4])
                        positions.append([x, y, z])
                        vx = float(L[5])
                        vy = float(L[6])
                        vz = float(L[7])
                        velocities.append([vx, vy, vz])
                        qw = float(L[8])
                        qx = float(L[9])
                        qy = float(L[10])
                        qz = float(L[11])
                        quaternions.append([qw, qx, qy, qz])
                        jx = float(L[12])
                        jy = float(L[13])
                        jz = float(L[14])
                        angVels.append([jx, jy, jz])
                        line = mdFile.readline()
                        if '</StuntDoubles>' in line:
                            break
                    break
        line = mdFile.readline()
        if not line: break
    
    mdFile.close()

def writeFile(outputFileName):
    outputFile = open(outputFileName, 'w')

    outputFile.write("<OOPSE version=4>\n");
    
    for metaline in metaData:
        outputFile.write(metaline)

    outputFile.write("  <Snapshot>\n")

    for frameline in frameData:
        outputFile.write(frameline)
        
    outputFile.write("    <StuntDoubles>\n")

    sdFormat = 'pvqj'
    for i in range(len(indices)):

        outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (indices[i], sdFormat, positions[i][0], positions[i][1], positions[i][2], velocities[i][0], velocities[i][1], velocities[i][2], quaternions[i][0], quaternions[i][1], quaternions[i][2], quaternions[i][3], angVels[i][0], angVels[i][1], angVels[i][2]))

    outputFile.write("    </StuntDoubles>\n")
    outputFile.write("  </Snapshot>\n")
    outputFile.write("</OOPSE>\n")
    outputFile.close()

def roundMe(x):
    if (x >= 0.0):
        return math.floor(x + 0.5)
    else:
        return math.ceil(x - 0.5)

def wrapVector(myVect):
    scaled = [0.0, 0.0, 0.0]
    for i in range(3):
        scaled[i] = myVect[i] * BoxInv[i]
        scaled[i] = scaled[i] - roundMe(scaled[i])
        myVect[i] = scaled[i] * Hmat[i][i]
    return myVect

def dot(L1, L2):
    myDot = 0.0
    for i in range(len(L1)):
        myDot = myDot + L1[i]*L2[i]
    return myDot

def normalize(L1):
    L2 = []
    myLength = math.sqrt(dot(L1, L1))
    for i in range(len(L1)):
        L2.append(L1[i] / myLength)
    return L2

def cross(L1, L2):
    # don't call this with anything other than length 3 lists please
    # or you'll be sorry
    L3 = [0.0, 0.0, 0.0]
    L3[0] = L1[1]*L2[2] - L1[2]*L2[1] 
    L3[1] = L1[2]*L2[0] - L1[0]*L2[2]
    L3[2] = L1[0]*L2[1] - L1[1]*L2[0]
    return L3

def f(x):
    return x

def findNeighbors():

    for i in range(len(indices)):
        neighbors.append(list())

    
    for i in range(len(indices)-1):
        iPos = positions[i]
        for j in range(i+1, len(indices)):
            jPos = positions[j]

            dpos = [jPos[0] - iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
            dpos = wrapVector(dpos)
            dist2 = dot(dpos,dpos)
            
            if (dist2 < 9.0):
                neighbors[i].append(j)
                neighbors[j].append(i)

            if (len(neighbors[i]) == 4):
                break

    surfaceCount = 0
    for i in range(len(indices)):
        if (len(neighbors[i]) == 3):
            neighbors[i].append(-1)
            surfaceCount = surfaceCount + 1

    print "surfaceCount = %d" % surfaceCount

def randomizeProtons():

    # do 10 proton bucket brigades:
    for i in range(10):
        origPoint = random.randint(0,len(indices))
        protonBucketBrigade(origPoint)

    # check to make sure everyone has a happy proton set:

    for j in range(len(indices)):
        if (len(DonatingTo[j]) != 2):
            # print "first round to fix molecule %d" % j
            protonBucketBrigade(j)

    # check to make sure everyone has a happy proton set:

    for j in range(len(indices)):
        if (len(DonatingTo[j]) != 2):
            print "second round to fix molecule %d" % j
            protonBucketBrigade(j)

    for j in range(len(indices)):
        if (len(DonatingTo[j]) != 2):
            print "unfilled proton donor list for molecule %d" % j


def protonBucketBrigade(origPoint):

    for i in range(len(indices)):
        DonatingTo.append(list())
        AcceptingFrom.append(list())
        OldDonor.append(list())

    # print "randomizing Quaternions"
    nIdeal = len(idealQuats)

    donor = origPoint
    # pick unreasonable start values (that don't match surface atom unreasonable values)
    acceptor = -10
    #oldDonor = -10
    # print 'origPoint = %d' % (origPoint)

    for i in range(50000):
    #while (acceptor != origPoint):
        myNeighbors = Set(neighbors[donor])
             
        # can't pick a proton choice from one of my current protons:
        badChoices = Set(DonatingTo[donor]).union(Set(OldDonor[acceptor]))
        # can't send a proton back to guy who sent one to me (no give-backs):
        #badChoices.add(Set(OldDonor[acceptor]))
        # can't send a proton to anyone who is already taking 2:
        for j in myNeighbors.difference(surfaceSet):
            if (len(AcceptingFrom[j]) == 2):
                badChoices.add(j)
                
        nDonors = len(DonatingTo[donor])

        if (nDonors <= 1):
            acceptor = random.choice(list( myNeighbors.difference(badChoices) ) )
            DonatingTo[donor].append(acceptor)
            if (acceptor != -1):
                AcceptingFrom[acceptor].append(donor)
                OldDonor[acceptor].append(donor)
        elif (nDonors == 2):
            acceptor = random.choice(DonatingTo[donor])
        else:
            print "Whoah!  How'd we get here?"

        #OldDonor = donor
        donor = acceptor
        if (acceptor == -1):
            # surface atoms all have a -1 neighbor, but that's OK.  A proton
            # is allowed to point out of the surface, but it does break the
            # proton chain letter
            # print "surface atom found, starting over from origPoint"
            donor = origPoint
            # break

def computeQuats():

    for i in range(len(indices)):
        DonatingTo.append(list())
        AcceptingFrom.append(list())
        OldDonor.append(list())
    # print "Computing Quaternions"
    ux = [0.0, 0.0, 0.0]
    uy = [0.0, 0.0, 0.0]
    uz = [0.0, 0.0, 0.0]
    RotMat = [ux, uy, uz]
    totalDipole = [0.0, 0.0, 0.0]
    for i in range(len(indices)):
        # print "doing quats for molecule %d" % i
        # print "Dlen = %d " % len(DonatingTo[i])
        # print DonatingTo[i]
        
        
        myPos = positions[i]

        acceptor1 = DonatingTo[i][0]
        if (acceptor1 == -1):
            tempVec = [0.0, 0.0, 0.0]
            for j in range(3):
                thisNeighbor = neighbors[i][j]
                npos = positions[thisNeighbor]
                npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
                npos1 = wrapVector(npos1)
                tempVec[0] = tempVec[0]  + npos1[0]
                tempVec[1] = tempVec[1]  + npos1[1]
                tempVec[2] = tempVec[2]  + npos1[2]                
            dpos1 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
            dpos1 = normalize(dpos1)
        else:
            a1pos = positions[acceptor1]
            dpos1 = [a1pos[0] - myPos[0], a1pos[1] - myPos[1], a1pos[2] - myPos[2]]
            dpos1 = wrapVector(dpos1)
            dpos1 = normalize(dpos1)
        
        acceptor2 = DonatingTo[i][1]
        if (acceptor2 == -1):
            tempVec = [0.0, 0.0, 0.0]
            for j in range(3):
                thisNeighbor = neighbors[i][j]
                npos = positions[thisNeighbor]
                npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
                npos1 = wrapVector(npos1)
                tempVec[0] = tempVec[0]  + npos1[0]
                tempVec[1] = tempVec[1]  + npos1[1]
                tempVec[2] = tempVec[2]  + npos1[2]                
            dpos2 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
            dpos2 = normalize(dpos2)
        else:
            a2pos = positions[acceptor2]
            dpos2 = [a2pos[0] - myPos[0], a2pos[1] - myPos[1], a2pos[2] - myPos[2]]
            dpos2 = wrapVector(dpos2)
            dpos2 = normalize(dpos2)


        for j in range(3):
            uz[j] = (dpos1[j] + dpos2[j])/2.0
        uz = normalize(uz)
        for j in range(3):
            uy[j] = dpos2[j] - dpos1[j]
        uy = normalize(uy)
        ux = cross(uy, uz)
        ux = normalize(ux)

        q = [0.0, 0.0, 0.0, 0.0]

        # RotMat to Quat code is out of OOPSE's SquareMatrix3.hpp code:

        RotMat[0] = ux
        RotMat[1] = uy
        RotMat[2] = uz

        t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
        
        if( t > 1e-6 ):
            s = 0.5 / math.sqrt( t )
            q[0] = 0.25 / s
            q[1] = (RotMat[1][2] - RotMat[2][1]) * s
            q[2] = (RotMat[2][0] - RotMat[0][2]) * s
            q[3] = (RotMat[0][1] - RotMat[1][0]) * s
        else:
            ad1 = RotMat[0][0]
            ad2 = RotMat[1][1]
            ad3 = RotMat[2][2]

            if( ad1 >= ad2 and ad1 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
                q[0] = (RotMat[1][2] - RotMat[2][1]) * s
                q[1] = 0.25 / s
                q[2] = (RotMat[0][1] + RotMat[1][0]) * s
                q[3] = (RotMat[0][2] + RotMat[2][0]) * s
            elif ( ad2 >= ad1 and ad2 >= ad3 ):
                s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
                q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
                q[1] = (RotMat[0][1] + RotMat[1][0]) * s
                q[2] = 0.25 / s
                q[3] = (RotMat[1][2] + RotMat[2][1]) * s
            else:
                s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
                q[0] = (RotMat[0][1] - RotMat[1][0]) * s
                q[1] = (RotMat[0][2] + RotMat[2][0]) * s
                q[2] = (RotMat[1][2] + RotMat[2][1]) * s
                q[3] = 0.25 / s

        quaternions[i] = q
        totalDipole = [totalDipole[0] + uz[0], totalDipole[1] + uz[1], totalDipole[2] + uz[2]]
    totalDipole = [-totalDipole[0], -totalDipole[1], -totalDipole[2]]
    print totalDipole
    Dipole = math.sqrt(dot(totalDipole, totalDipole))
    print 'Total Dipole Moment = %d' % Dipole
    if (Dipole > 20 or Dipole < -20):
        print "Bad Dipole, starting over"
        for i in range(len(indices)):
            del OldDonor[:]
            del AcceptingFrom[:]
            del DonatingTo[:]
            # del badChoices[:]
        randomizeProtons()
        computeQuats()
    else:
        print "All Done!"
        

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hm:o:", ["help", "meta-data=", "output-file="]) 
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()                  
        elif opt in ("-m", "--meta-data"): 
            mdFileName = arg
            global _haveMDFileName
            _haveMDFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            global _haveOutputFileName
            _haveOutputFileName = 1
    if (_haveMDFileName != 1):
        usage() 
        print "No meta-data file was specified"
        sys.exit()
    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()
    readFile(mdFileName)
    # analyzeQuats()
    findNeighbors()
    randomizeProtons()
    computeQuats()
    writeFile(outputFileName)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1189
Committed:	Mon Nov 26 19:23:36 2007 UTC (17 years, 11 months ago) by cpuglis
Original Path:	*trunk/src/applications/utilities/waterRotator*
File size:	16041 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	cpuglis	1189	#!/usr/bin/env python
2			"""Water Quaternion Sampler
3
4			Samples water orientations from a list of known good
5			orientations
6
7			Usage: waterRotator
8
9			Options:
10			-h, --help show this help
11			-m, --meta-data=... use specified meta-data (.md) file
12			-o, --output-file=... use specified output (.md) file
13
14
15			Example:
16			waterRotator -m basal.md -o basal.new.md
17
18			"""
19
20			__author__ = "Dan Gezelter (gezelter@nd.edu)"
21			__version__ = "$Revision: 1.1 $"
22			__date__ = "$Date: 2007-11-26 19:23:36 $"
23			__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
24			__license__ = "OOPSE"
25
26			import sys
27			import getopt
28			import string
29			import math
30			import random
31			from sets import *
32			#from Numeric import *
33
34			_haveMDFileName = 0
35			_haveOutputFileName = 0
36
37			metaData = []
38			frameData = []
39			positions = []
40			velocities = []
41			quaternions = []
42			angVels = []
43			indices = []
44			neighbors = []
45			DonatingTo = []
46			AcceptingFrom = []
47			OldDonor = []
48			Hmat = []
49			BoxInv = []
50			#Hmat = zeros([3,3],Float)
51			#BoxInv = zeros([3],Float)
52			H1vects = []
53			H2vects = []
54			DipoleVects = []
55			allAcceptors = []
56			availableneighbs = []
57			surfaceSet = Set([-1])
58
59
60
61			def usage():
62			print __doc__
63
64
65			def readFile(mdFileName):
66			mdFile = open(mdFileName, 'r')
67			# Find OOPSE version info first
68			line = mdFile.readline()
69			while 1:
70			if '<OOPSE version=' in line:
71			OOPSEversion = line
72			break
73			line = mdFile.readline()
74
75			# Rewind file and find start of MetaData block
76
77			mdFile.seek(0)
78			line = mdFile.readline()
79			print "reading MetaData"
80			while 1:
81			if '<MetaData>' in line:
82			while 2:
83			metaData.append(line)
84			line = mdFile.readline()
85			if '</MetaData>' in line:
86			metaData.append(line)
87			break
88			break
89			line = mdFile.readline()
90
91			mdFile.seek(0)
92			print "reading Snapshot"
93			line = mdFile.readline()
94			while 1:
95			if '<Snapshot>' in line:
96			line = mdFile.readline()
97			while 1:
98			print "reading FrameData"
99			if '<FrameData>' in line:
100			while 2:
101			frameData.append(line)
102			if 'Hmat:' in line:
103			L = line.split()
104			Hxx = float(L[2].strip(','))
105			Hxy = float(L[3].strip(','))
106			Hxz = float(L[4].strip(','))
107			Hyx = float(L[7].strip(','))
108			Hyy = float(L[8].strip(','))
109			Hyz = float(L[9].strip(','))
110			Hzx = float(L[12].strip(','))
111			Hzy = float(L[13].strip(','))
112			Hzz = float(L[14].strip(','))
113			Hmat.append([Hxx, Hxy, Hxz])
114			Hmat.append([Hyx, Hyy, Hyz])
115			Hmat.append([Hzx, Hzy, Hzz])
116			print Hmat
117			BoxInv.append(1.0/Hxx)
118			BoxInv.append(1.0/Hyy)
119			BoxInv.append(1.0/Hzz)
120			print BoxInv
121			line = mdFile.readline()
122			if '</FrameData>' in line:
123			frameData.append(line)
124			break
125			break
126
127			line = mdFile.readline()
128			while 1:
129			if '<StuntDoubles>' in line:
130			line = mdFile.readline()
131			while 2:
132			L = line.split()
133			myIndex = int(L[0])
134			indices.append(myIndex)
135			pvqj = L[1]
136			x = float(L[2])
137			y = float(L[3])
138			z = float(L[4])
139			positions.append([x, y, z])
140			vx = float(L[5])
141			vy = float(L[6])
142			vz = float(L[7])
143			velocities.append([vx, vy, vz])
144			qw = float(L[8])
145			qx = float(L[9])
146			qy = float(L[10])
147			qz = float(L[11])
148			quaternions.append([qw, qx, qy, qz])
149			jx = float(L[12])
150			jy = float(L[13])
151			jz = float(L[14])
152			angVels.append([jx, jy, jz])
153			line = mdFile.readline()
154			if '</StuntDoubles>' in line:
155			break
156			break
157			line = mdFile.readline()
158			if not line: break
159
160			mdFile.close()
161
162			def writeFile(outputFileName):
163			outputFile = open(outputFileName, 'w')
164
165			outputFile.write("<OOPSE version=4>\n");
166
167			for metaline in metaData:
168			outputFile.write(metaline)
169
170			outputFile.write(" <Snapshot>\n")
171
172			for frameline in frameData:
173			outputFile.write(frameline)
174
175			outputFile.write(" <StuntDoubles>\n")
176
177			sdFormat = 'pvqj'
178			for i in range(len(indices)):
179
180			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (indices[i], sdFormat, positions[i][0], positions[i][1], positions[i][2], velocities[i][0], velocities[i][1], velocities[i][2], quaternions[i][0], quaternions[i][1], quaternions[i][2], quaternions[i][3], angVels[i][0], angVels[i][1], angVels[i][2]))
181
182			outputFile.write(" </StuntDoubles>\n")
183			outputFile.write(" </Snapshot>\n")
184			outputFile.write("</OOPSE>\n")
185			outputFile.close()
186
187			def roundMe(x):
188			if (x >= 0.0):
189			return math.floor(x + 0.5)
190			else:
191			return math.ceil(x - 0.5)
192
193			def wrapVector(myVect):
194			scaled = [0.0, 0.0, 0.0]
195			for i in range(3):
196			scaled[i] = myVect[i] * BoxInv[i]
197			scaled[i] = scaled[i] - roundMe(scaled[i])
198			myVect[i] = scaled[i] * Hmat[i][i]
199			return myVect
200
201			def dot(L1, L2):
202			myDot = 0.0
203			for i in range(len(L1)):
204			myDot = myDot + L1[i]*L2[i]
205			return myDot
206
207			def normalize(L1):
208			L2 = []
209			myLength = math.sqrt(dot(L1, L1))
210			for i in range(len(L1)):
211			L2.append(L1[i] / myLength)
212			return L2
213
214			def cross(L1, L2):
215			# don't call this with anything other than length 3 lists please
216			# or you'll be sorry
217			L3 = [0.0, 0.0, 0.0]
218			L3[0] = L1[1]L2[2] - L1[2]L2[1]
219			L3[1] = L1[2]L2[0] - L1[0]L2[2]
220			L3[2] = L1[0]L2[1] - L1[1]L2[0]
221			return L3
222
223			def f(x):
224			return x
225
226			def findNeighbors():
227
228			for i in range(len(indices)):
229			neighbors.append(list())
230
231
232			for i in range(len(indices)-1):
233			iPos = positions[i]
234			for j in range(i+1, len(indices)):
235			jPos = positions[j]
236
237			dpos = [jPos[0] - iPos[0], jPos[1]-iPos[1], jPos[2]-iPos[2]]
238			dpos = wrapVector(dpos)
239			dist2 = dot(dpos,dpos)
240
241			if (dist2 < 9.0):
242			neighbors[i].append(j)
243			neighbors[j].append(i)
244
245			if (len(neighbors[i]) == 4):
246			break
247
248			surfaceCount = 0
249			for i in range(len(indices)):
250			if (len(neighbors[i]) == 3):
251			neighbors[i].append(-1)
252			surfaceCount = surfaceCount + 1
253
254			print "surfaceCount = %d" % surfaceCount
255
256			def randomizeProtons():
257
258			# do 10 proton bucket brigades:
259			for i in range(10):
260			origPoint = random.randint(0,len(indices))
261			protonBucketBrigade(origPoint)
262
263			# check to make sure everyone has a happy proton set:
264
265			for j in range(len(indices)):
266			if (len(DonatingTo[j]) != 2):
267			# print "first round to fix molecule %d" % j
268			protonBucketBrigade(j)
269
270			# check to make sure everyone has a happy proton set:
271
272			for j in range(len(indices)):
273			if (len(DonatingTo[j]) != 2):
274			print "second round to fix molecule %d" % j
275			protonBucketBrigade(j)
276
277			for j in range(len(indices)):
278			if (len(DonatingTo[j]) != 2):
279			print "unfilled proton donor list for molecule %d" % j
280
281
282
283			def protonBucketBrigade(origPoint):
284
285			for i in range(len(indices)):
286			DonatingTo.append(list())
287			AcceptingFrom.append(list())
288			OldDonor.append(list())
289
290			# print "randomizing Quaternions"
291			nIdeal = len(idealQuats)
292
293			donor = origPoint
294			# pick unreasonable start values (that don't match surface atom unreasonable values)
295			acceptor = -10
296			#oldDonor = -10
297			# print 'origPoint = %d' % (origPoint)
298
299			for i in range(50000):
300			#while (acceptor != origPoint):
301			myNeighbors = Set(neighbors[donor])
302
303			# can't pick a proton choice from one of my current protons:
304			badChoices = Set(DonatingTo[donor]).union(Set(OldDonor[acceptor]))
305			# can't send a proton back to guy who sent one to me (no give-backs):
306			#badChoices.add(Set(OldDonor[acceptor]))
307			# can't send a proton to anyone who is already taking 2:
308			for j in myNeighbors.difference(surfaceSet):
309			if (len(AcceptingFrom[j]) == 2):
310			badChoices.add(j)
311
312			nDonors = len(DonatingTo[donor])
313
314			if (nDonors <= 1):
315			acceptor = random.choice(list( myNeighbors.difference(badChoices) ) )
316			DonatingTo[donor].append(acceptor)
317			if (acceptor != -1):
318			AcceptingFrom[acceptor].append(donor)
319			OldDonor[acceptor].append(donor)
320			elif (nDonors == 2):
321			acceptor = random.choice(DonatingTo[donor])
322			else:
323			print "Whoah! How'd we get here?"
324
325			#OldDonor = donor
326			donor = acceptor
327			if (acceptor == -1):
328			# surface atoms all have a -1 neighbor, but that's OK. A proton
329			# is allowed to point out of the surface, but it does break the
330			# proton chain letter
331			# print "surface atom found, starting over from origPoint"
332			donor = origPoint
333			# break
334
335			def computeQuats():
336
337			for i in range(len(indices)):
338			DonatingTo.append(list())
339			AcceptingFrom.append(list())
340			OldDonor.append(list())
341			# print "Computing Quaternions"
342			ux = [0.0, 0.0, 0.0]
343			uy = [0.0, 0.0, 0.0]
344			uz = [0.0, 0.0, 0.0]
345			RotMat = [ux, uy, uz]
346			totalDipole = [0.0, 0.0, 0.0]
347			for i in range(len(indices)):
348			# print "doing quats for molecule %d" % i
349			# print "Dlen = %d " % len(DonatingTo[i])
350			# print DonatingTo[i]
351
352
353			myPos = positions[i]
354
355			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
374			acceptor2 = DonatingTo[i][1]
375			if (acceptor2 == -1):
376			tempVec = [0.0, 0.0, 0.0]
377			for j in range(3):
378			thisNeighbor = neighbors[i][j]
379			npos = positions[thisNeighbor]
380			npos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
381			npos1 = wrapVector(npos1)
382			tempVec[0] = tempVec[0] + npos1[0]
383			tempVec[1] = tempVec[1] + npos1[1]
384			tempVec[2] = tempVec[2] + npos1[2]
385			dpos2 = [-tempVec[0]/3.0, -tempVec[1]/3.0, -tempVec[2]/3.0]
386			dpos2 = normalize(dpos2)
387			else:
388			a2pos = positions[acceptor2]
389			dpos2 = [a2pos[0] - myPos[0], a2pos[1] - myPos[1], a2pos[2] - myPos[2]]
390			dpos2 = wrapVector(dpos2)
391			dpos2 = normalize(dpos2)
392
393
394
395
396
397			for j in range(3):
398			uz[j] = (dpos1[j] + dpos2[j])/2.0
399			uz = normalize(uz)
400			for j in range(3):
401			uy[j] = dpos2[j] - dpos1[j]
402			uy = normalize(uy)
403			ux = cross(uy, uz)
404			ux = normalize(ux)
405
406			q = [0.0, 0.0, 0.0, 0.0]
407
408			# RotMat to Quat code is out of OOPSE's SquareMatrix3.hpp code:
409
410			RotMat[0] = ux
411			RotMat[1] = uy
412			RotMat[2] = uz
413
414			t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
415
416			if( t > 1e-6 ):
417			s = 0.5 / math.sqrt( t )
418			q[0] = 0.25 / s
419			q[1] = (RotMat[1][2] - RotMat[2][1]) * s
420			q[2] = (RotMat[2][0] - RotMat[0][2]) * s
421			q[3] = (RotMat[0][1] - RotMat[1][0]) * s
422			else:
423			ad1 = RotMat[0][0]
424			ad2 = RotMat[1][1]
425			ad3 = RotMat[2][2]
426
427			if( ad1 >= ad2 and ad1 >= ad3 ):
428			s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
429			q[0] = (RotMat[1][2] - RotMat[2][1]) * s
430			q[1] = 0.25 / s
431			q[2] = (RotMat[0][1] + RotMat[1][0]) * s
432			q[3] = (RotMat[0][2] + RotMat[2][0]) * s
433			elif ( ad2 >= ad1 and ad2 >= ad3 ):
434			s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
435			q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
436			q[1] = (RotMat[0][1] + RotMat[1][0]) * s
437			q[2] = 0.25 / s
438			q[3] = (RotMat[1][2] + RotMat[2][1]) * s
439			else:
440			s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
441			q[0] = (RotMat[0][1] - RotMat[1][0]) * s
442			q[1] = (RotMat[0][2] + RotMat[2][0]) * s
443			q[2] = (RotMat[1][2] + RotMat[2][1]) * s
444			q[3] = 0.25 / s
445
446			quaternions[i] = q
447			totalDipole = [totalDipole[0] + uz[0], totalDipole[1] + uz[1], totalDipole[2] + uz[2]]
448			totalDipole = [-totalDipole[0], -totalDipole[1], -totalDipole[2]]
449			print totalDipole
450			Dipole = math.sqrt(dot(totalDipole, totalDipole))
451			print 'Total Dipole Moment = %d' % Dipole
452			if (Dipole > 20 or Dipole < -20):
453			print "Bad Dipole, starting over"
454			for i in range(len(indices)):
455			del OldDonor[:]
456			del AcceptingFrom[:]
457			del DonatingTo[:]
458			# del badChoices[:]
459			randomizeProtons()
460			computeQuats()
461			else:
462			print "All Done!"
463
464
465
466			def main(argv):
467			try:
468			opts, args = getopt.getopt(argv, "hm:o:", ["help", "meta-data=", "output-file="])
469			except getopt.GetoptError:
470			usage()
471			sys.exit(2)
472			for opt, arg in opts:
473			if opt in ("-h", "--help"):
474			usage()
475			sys.exit()
476			elif opt in ("-m", "--meta-data"):
477			mdFileName = arg
478			global _haveMDFileName
479			_haveMDFileName = 1
480			elif opt in ("-o", "--output-file"):
481			outputFileName = arg
482			global _haveOutputFileName
483			_haveOutputFileName = 1
484			if (_haveMDFileName != 1):
485			usage()
486			print "No meta-data file was specified"
487			sys.exit()
488			if (_haveOutputFileName != 1):
489			usage()
490			print "No output file was specified"
491			sys.exit()
492			readFile(mdFileName)
493			# analyzeQuats()
494			findNeighbors()
495			randomizeProtons()
496			computeQuats()
497			writeFile(outputFileName)
498
499			if __name__ == "__main__":
500			if len(sys.argv) == 1:
501			usage()
502			sys.exit()
503			main(sys.argv[1:])