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

#!@PYTHON_EXECUTABLE@
"""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$"
__date__ = "$Date$"
__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import random

_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 OpenMD version info first
    line = mdFile.readline()
    while 1:
        if '<OOPSE version=' in line or '<OpenMD version=' in line:
            OpenMDversion = 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("<OpenMD version=1>\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("</OpenMD>\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]) < 4):
            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())

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

    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):

            if (len(myNeighbors.difference(badChoices)) != 0):
                acceptor = random.choice(list(myNeighbors.difference(badChoices)))
            else:
                acceptor = -1

            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?"

        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]
        acceptor2 = DonatingTo[i][1]
        if (acceptor1 == -1 and acceptor2 == -1):
            donor1 = AcceptingFrom[i][0]            
            donor2 = AcceptingFrom[i][1]
            npos = positions[donor1]
            apos1 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
            apos1 = wrapVector(apos1)
            npos = positions[donor2]
            apos2 = [npos[0] - myPos[0], npos[1] - myPos[1], npos[2] - myPos[2]]
            apos2 = wrapVector(apos2)
            tempX = [0.0, 0.0, 0.0]
            tempY = [0.0, 0.0, 0.0]
            tempZ = [0.0, 0.0, 0.0]
            tempZ[0] = 0.5*(apos1[0] + apos2[0])
            tempZ[1] = 0.5*(apos1[1] + apos2[1])
            tempZ[2] = 0.5*(apos1[2] + apos2[2])
            tempX[0] = (apos2[0] - apos1[0])
            tempX[1] = (apos2[1] - apos1[1])
            tempX[2] = (apos2[2] - apos1[2])
            tempZ = normalize(tempZ)
            tempX = normalize(tempX)
            tempY = cross(tempX, tempZ)
            dpos1[0] = tempZ[0] - 0.5*tempY[0]
            dpos1[1] = tempZ[1] - 0.5*tempY[1]
            dpos1[2] = tempZ[2] - 0.5*tempY[2]
            dpos2[0] = tempZ[0] + 0.5*tempY[0]
            dpos2[1] = tempZ[1] + 0.5*tempY[1]
            dpos2[2] = tempZ[2] + 0.5*tempY[2]
            
        else:
            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)
        

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