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

#!/usr/bin/env python
"""Water Replacer

Finds atomistic waters in an xyz file and generates a meta-data file
with center of mass and orientational coordinates for rigid body
waters.

Usage: waterReplacer

Options:
  -h, --help              show this help
  -x,                     use the specified input (.xyz) file
  -o, --output-file=...   use specified output (.md) file


Example:
   waterReplacer -x basal.xyz -o basal.md

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1646 $"
__date__ = "$Date: 2011-09-26 09:30:00 -0400 (Mon, 26 Sep 2011) $"
__copyright__ = "Copyright (c) 2011 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math
import random
import numpy

_haveXYZFileName = 0
_haveOutputFileName = 0

atypes = []
positions = []
metaData = []
frameData = []
WaterPos = []
WaterQuats = []
indices = []
Hmat = []
BoxInv = []
H = []
Eliminate = []

#Hmat = zeros([3,3],Float)
#BoxInv = zeros([3],Float)

def usage():
    print __doc__

def readFile(XYZFileName):
    print "reading XYZ file"

    XYZFile = open(XYZFileName, 'r')        
    # Find number of atoms first
    line = XYZFile.readline()
    L = line.split()
    nAtoms = int(L[0])
    # skip comment line
    line = XYZFile.readline()
    for i in range(nAtoms):
        line = XYZFile.readline()
        L = line.split()
        myIndex = i
        indices.append(myIndex)
        atomType = L[0]
        atypes.append(atomType)
        x = float(L[1])
        y = float(L[2])
        z = float(L[3])
        positions.append([x, y, z])
    XYZFile.close()

def findWaters():
    print "finding water molecules"
    # simpler since we only have to find H atoms within a few
    # angstroms of each water:

    hCovRad = 0.32
    oCovRad = 0.73
    covTol = 0.45
    OHbond = hCovRad + oCovRad + covTol
    Hmass = 1.0079
    Omass = 15.9994
    # initialize H array to an error condition
    H.append(-1)
    H.append(-1)
    
    for i in range(len(indices)):
        if (atypes[i] == "O"):
            nH = 0
            H[0] = -1
            H[1] = -1
            COM = [0.0, 0.0, 0.0]
            opos = positions[i]
            for j in range(len(indices)):
                if (atypes[j] == "H"):
                    hpos = positions[j]
                    dx = opos[0] - hpos[0]
                    dy = opos[1] - hpos[1]
                    dz = opos[2] - hpos[2]
                    dist = math.sqrt(dx*dx + dy*dy + dz*dz)
                    if (dist < OHbond):
                        if (nH >= 2):
                            print "oxygen %d had too many hydrogens" % (i)
                            sys.exit(1)
                        H[nH] = j
                        nH = nH + 1
            if (nH != 2):
                print "oxygen %d had %d hydrogens, skipping" % (i, nH)
            if (nH == 2):
                Xcom = Omass * opos[0] + Hmass*(positions[H[0]][0] + positions[H[1]][0])
                Ycom = Omass * opos[1] + Hmass*(positions[H[0]][1] + positions[H[1]][1])
                Zcom = Omass * opos[2] + Hmass*(positions[H[0]][2] + positions[H[1]][2])
                
                totalMass = Omass + 2.0*Hmass
                Xcom = Xcom / totalMass
                Ycom = Ycom / totalMass
                Zcom = Zcom / totalMass
                COM = [Xcom, Ycom, Zcom]
                WaterPos.append(COM)
                bisector = [0.0, 0.0, 0.0]
                ux = [0.0, 0.0, 0.0]
                uy = [0.0, 0.0, 0.0]
                uz = [0.0, 0.0, 0.0]
                RotMat = numpy.zeros((3,3), numpy.float)

                for j in range(3):
                    bisector[j] = 0.5*(positions[H[0]][j] + positions[H[1]][j])
                    uz[j] = bisector[j] - opos[j]
                    uy[j] = positions[H[0]][j] - positions[H[1]][j]
                    
                uz = normalize(uz)
                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

                WaterQuats.append(q)
                Eliminate.append(i)
                Eliminate.append(H[0])
                Eliminate.append(H[1])


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'

    index = 0
    for i in range(len(WaterPos)):
        outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (index, sdFormat, WaterPos[i][0], WaterPos[i][1], WaterPos[i][2], 0.0, 0.0, 0.0, WaterQuats[i][0], WaterQuats[i][1], WaterQuats[i][2], WaterQuats[i][3], 0.0, 0.0, 0.0))
        index = index + 1


    sdFormat = 'pv'
    for i in range(len(indices)):
        if  i not in Eliminate:
            outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e \n" % (index, sdFormat, positions[i][0], positions[i][1], positions[i][2], 0.0, 0.0, 0.0))

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

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 main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz-file=", "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 ("-x", "--xyz-file"): 
            xyzFileName = arg
            global _haveXYZFileName
            _haveXYZFileName = 1
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            global _haveOutputFileName
            _haveOutputFileName = 1
    if (_haveXYZFileName != 1):
        usage() 
        print "No XYZ file was specified"
        sys.exit()
    if (_haveOutputFileName != 1):
        usage()
        print "No output file was specified"
        sys.exit()
    readFile(xyzFileName)
    findWaters()
    writeFile(outputFileName)

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1663
Committed:	Tue Nov 15 16:57:53 2011 UTC (13 years, 11 months ago) by gezelter
File size:	8870 byte(s)
Log Message:	Added a waterReplacer script. A bit rough around the edges, but seems to do the right thing for a slab of ice.
#	User	Rev	Content
1	gezelter	1663	#!/usr/bin/env python
2			"""Water Replacer
3
4			Finds atomistic waters in an xyz file and generates a meta-data file
5			with center of mass and orientational coordinates for rigid body
6			waters.
7
8			Usage: waterReplacer
9
10			Options:
11			-h, --help show this help
12			-x, use the specified input (.xyz) file
13			-o, --output-file=... use specified output (.md) file
14
15
16			Example:
17			waterReplacer -x basal.xyz -o basal.md
18
19			"""
20
21			__author__ = "Dan Gezelter (gezelter@nd.edu)"
22			__version__ = "$Revision: 1646 $"
23			__date__ = "$Date: 2011-09-26 09:30:00 -0400 (Mon, 26 Sep 2011) $"
24			__copyright__ = "Copyright (c) 2011 by the University of Notre Dame"
25			__license__ = "OpenMD"
26
27			import sys
28			import getopt
29			import string
30			import math
31			import random
32			import numpy
33
34			_haveXYZFileName = 0
35			_haveOutputFileName = 0
36
37			atypes = []
38			positions = []
39			metaData = []
40			frameData = []
41			WaterPos = []
42			WaterQuats = []
43			indices = []
44			Hmat = []
45			BoxInv = []
46			H = []
47			Eliminate = []
48
49			#Hmat = zeros([3,3],Float)
50			#BoxInv = zeros([3],Float)
51
52			def usage():
53			print __doc__
54
55			def readFile(XYZFileName):
56			print "reading XYZ file"
57
58			XYZFile = open(XYZFileName, 'r')
59			# Find number of atoms first
60			line = XYZFile.readline()
61			L = line.split()
62			nAtoms = int(L[0])
63			# skip comment line
64			line = XYZFile.readline()
65			for i in range(nAtoms):
66			line = XYZFile.readline()
67			L = line.split()
68			myIndex = i
69			indices.append(myIndex)
70			atomType = L[0]
71			atypes.append(atomType)
72			x = float(L[1])
73			y = float(L[2])
74			z = float(L[3])
75			positions.append([x, y, z])
76			XYZFile.close()
77
78			def findWaters():
79			print "finding water molecules"
80			# simpler since we only have to find H atoms within a few
81			# angstroms of each water:
82
83			hCovRad = 0.32
84			oCovRad = 0.73
85			covTol = 0.45
86			OHbond = hCovRad + oCovRad + covTol
87			Hmass = 1.0079
88			Omass = 15.9994
89			# initialize H array to an error condition
90			H.append(-1)
91			H.append(-1)
92
93			for i in range(len(indices)):
94			if (atypes[i] == "O"):
95			nH = 0
96			H[0] = -1
97			H[1] = -1
98			COM = [0.0, 0.0, 0.0]
99			opos = positions[i]
100			for j in range(len(indices)):
101			if (atypes[j] == "H"):
102			hpos = positions[j]
103			dx = opos[0] - hpos[0]
104			dy = opos[1] - hpos[1]
105			dz = opos[2] - hpos[2]
106			dist = math.sqrt(dxdx + dydy + dz*dz)
107			if (dist < OHbond):
108			if (nH >= 2):
109			print "oxygen %d had too many hydrogens" % (i)
110			sys.exit(1)
111			H[nH] = j
112			nH = nH + 1
113			if (nH != 2):
114			print "oxygen %d had %d hydrogens, skipping" % (i, nH)
115			if (nH == 2):
116			Xcom = Omass * opos[0] + Hmass*(positions[H[0]][0] + positions[H[1]][0])
117			Ycom = Omass * opos[1] + Hmass*(positions[H[0]][1] + positions[H[1]][1])
118			Zcom = Omass * opos[2] + Hmass*(positions[H[0]][2] + positions[H[1]][2])
119
120			totalMass = Omass + 2.0*Hmass
121			Xcom = Xcom / totalMass
122			Ycom = Ycom / totalMass
123			Zcom = Zcom / totalMass
124			COM = [Xcom, Ycom, Zcom]
125			WaterPos.append(COM)
126			bisector = [0.0, 0.0, 0.0]
127			ux = [0.0, 0.0, 0.0]
128			uy = [0.0, 0.0, 0.0]
129			uz = [0.0, 0.0, 0.0]
130			RotMat = numpy.zeros((3,3), numpy.float)
131
132			for j in range(3):
133			bisector[j] = 0.5*(positions[H[0]][j] + positions[H[1]][j])
134			uz[j] = bisector[j] - opos[j]
135			uy[j] = positions[H[0]][j] - positions[H[1]][j]
136
137			uz = normalize(uz)
138			uy = normalize(uy)
139			ux = cross(uy, uz)
140			ux = normalize(ux)
141
142			q = [0.0, 0.0, 0.0, 0.0]
143
144			# RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
145
146			RotMat[0] = ux
147			RotMat[1] = uy
148			RotMat[2] = uz
149
150			t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
151
152			if( t > 1e-6 ):
153			s = 0.5 / math.sqrt( t )
154			q[0] = 0.25 / s
155			q[1] = (RotMat[1][2] - RotMat[2][1]) * s
156			q[2] = (RotMat[2][0] - RotMat[0][2]) * s
157			q[3] = (RotMat[0][1] - RotMat[1][0]) * s
158			else:
159			ad1 = RotMat[0][0]
160			ad2 = RotMat[1][1]
161			ad3 = RotMat[2][2]
162
163			if( ad1 >= ad2 and ad1 >= ad3 ):
164			s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
165			q[0] = (RotMat[1][2] - RotMat[2][1]) * s
166			q[1] = 0.25 / s
167			q[2] = (RotMat[0][1] + RotMat[1][0]) * s
168			q[3] = (RotMat[0][2] + RotMat[2][0]) * s
169			elif ( ad2 >= ad1 and ad2 >= ad3 ):
170			s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
171			q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
172			q[1] = (RotMat[0][1] + RotMat[1][0]) * s
173			q[2] = 0.25 / s
174			q[3] = (RotMat[1][2] + RotMat[2][1]) * s
175			else:
176			s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
177			q[0] = (RotMat[0][1] - RotMat[1][0]) * s
178			q[1] = (RotMat[0][2] + RotMat[2][0]) * s
179			q[2] = (RotMat[1][2] + RotMat[2][1]) * s
180			q[3] = 0.25 / s
181
182			WaterQuats.append(q)
183			Eliminate.append(i)
184			Eliminate.append(H[0])
185			Eliminate.append(H[1])
186
187
188			def writeFile(outputFileName):
189			outputFile = open(outputFileName, 'w')
190
191			outputFile.write("<OpenMD version=1>\n");
192
193			for metaline in metaData:
194			outputFile.write(metaline)
195
196			outputFile.write(" <Snapshot>\n")
197
198			for frameline in frameData:
199			outputFile.write(frameline)
200
201			outputFile.write(" <StuntDoubles>\n")
202
203			sdFormat = 'pvqj'
204
205			index = 0
206			for i in range(len(WaterPos)):
207			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e %13e %13e %13e %13e %13e %13e %13e\n" % (index, sdFormat, WaterPos[i][0], WaterPos[i][1], WaterPos[i][2], 0.0, 0.0, 0.0, WaterQuats[i][0], WaterQuats[i][1], WaterQuats[i][2], WaterQuats[i][3], 0.0, 0.0, 0.0))
208			index = index + 1
209
210
211			sdFormat = 'pv'
212			for i in range(len(indices)):
213			if i not in Eliminate:
214			outputFile.write("%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e \n" % (index, sdFormat, positions[i][0], positions[i][1], positions[i][2], 0.0, 0.0, 0.0))
215
216			outputFile.write(" </StuntDoubles>\n")
217			outputFile.write(" </Snapshot>\n")
218			outputFile.write("</OpenMD>\n")
219			outputFile.close()
220
221			def dot(L1, L2):
222			myDot = 0.0
223			for i in range(len(L1)):
224			myDot = myDot + L1[i]*L2[i]
225			return myDot
226
227			def normalize(L1):
228			L2 = []
229			myLength = math.sqrt(dot(L1, L1))
230			for i in range(len(L1)):
231			L2.append(L1[i] / myLength)
232			return L2
233
234			def cross(L1, L2):
235			# don't call this with anything other than length 3 lists please
236			# or you'll be sorry
237			L3 = [0.0, 0.0, 0.0]
238			L3[0] = L1[1]L2[2] - L1[2]L2[1]
239			L3[1] = L1[2]L2[0] - L1[0]L2[2]
240			L3[2] = L1[0]L2[1] - L1[1]L2[0]
241			return L3
242
243			def main(argv):
244			try:
245			opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz-file=", "output-file="])
246			except getopt.GetoptError:
247			usage()
248			sys.exit(2)
249			for opt, arg in opts:
250			if opt in ("-h", "--help"):
251			usage()
252			sys.exit()
253			elif opt in ("-x", "--xyz-file"):
254			xyzFileName = arg
255			global _haveXYZFileName
256			_haveXYZFileName = 1
257			elif opt in ("-o", "--output-file"):
258			outputFileName = arg
259			global _haveOutputFileName
260			_haveOutputFileName = 1
261			if (_haveXYZFileName != 1):
262			usage()
263			print "No XYZ file was specified"
264			sys.exit()
265			if (_haveOutputFileName != 1):
266			usage()
267			print "No output file was specified"
268			sys.exit()
269			readFile(xyzFileName)
270			findWaters()
271			writeFile(outputFileName)
272
273			if __name__ == "__main__":
274			if len(sys.argv) == 1:
275			usage()
276			sys.exit()
277			main(sys.argv[1:])