[ViewVC] Contents of: OpenMD/trunk/src/applications/utilities/waterReplacer

#!@PYTHON_EXECUTABLE@
"""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))
            index = index + 1

    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:	1879
Committed:	Sun Jun 16 15:15:42 2013 UTC (11 years, 10 months ago) by gezelter
File size:	8900 byte(s)
Log Message:	MERGE OpenMD development 1783:1878 into trunk
#	Content
1	#!@PYTHON_EXECUTABLE@
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	index = index + 1
216
217	outputFile.write(" </StuntDoubles>\n")
218	outputFile.write(" </Snapshot>\n")
219	outputFile.write("</OpenMD>\n")
220	outputFile.close()
221
222	def dot(L1, L2):
223	myDot = 0.0
224	for i in range(len(L1)):
225	myDot = myDot + L1[i]*L2[i]
226	return myDot
227
228	def normalize(L1):
229	L2 = []
230	myLength = math.sqrt(dot(L1, L1))
231	for i in range(len(L1)):
232	L2.append(L1[i] / myLength)
233	return L2
234
235	def cross(L1, L2):
236	# don't call this with anything other than length 3 lists please
237	# or you'll be sorry
238	L3 = [0.0, 0.0, 0.0]
239	L3[0] = L1[1]L2[2] - L1[2]L2[1]
240	L3[1] = L1[2]L2[0] - L1[0]L2[2]
241	L3[2] = L1[0]L2[1] - L1[1]L2[0]
242	return L3
243
244	def main(argv):
245	try:
246	opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz-file=", "output-file="])
247	except getopt.GetoptError:
248	usage()
249	sys.exit(2)
250	for opt, arg in opts:
251	if opt in ("-h", "--help"):
252	usage()
253	sys.exit()
254	elif opt in ("-x", "--xyz-file"):
255	xyzFileName = arg
256	global _haveXYZFileName
257	_haveXYZFileName = 1
258	elif opt in ("-o", "--output-file"):
259	outputFileName = arg
260	global _haveOutputFileName
261	_haveOutputFileName = 1
262	if (_haveXYZFileName != 1):
263	usage()
264	print "No XYZ file was specified"
265	sys.exit()
266	if (_haveOutputFileName != 1):
267	usage()
268	print "No output file was specified"
269	sys.exit()
270	readFile(xyzFileName)
271	findWaters()
272	writeFile(outputFileName)
273
274	if __name__ == "__main__":
275	if len(sys.argv) == 1:
276	usage()
277	sys.exit()
278	main(sys.argv[1:])