[ViewVC] Annotation of: OpenMD/trunk/src/applications/utilities/principalAxisCalculator

#!/usr/bin/env python
"""principalAxisCalculator

Opens an XYZ file and computes the moments of inertia and principal axes
for the structure in the XYZ file.  Optionally rotate the structure so that
the long axis (that with the smallest eigenvalue) is pointing along the
z-axis.

Usage: principalAxisCalculator

Options:
  -h,  --help              show this help
  -x,  --xyz=...           use specified XYZ (.xyz) file for the structure
  -o,  --out=...           rotate the structure so that the smallest eigenvalue
                           of the rotation matrix points along the z-axis.

Example:
   principalAxisCalculator -x junk.xyz -o rot.xyz 

"""

__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
from sets import *
import numpy


_haveXYZFileName = 0
_haveOutFileName = 0

positions = []
indices = []
atypes = []
Hmass = 1.0079
Cmass = 12.011
Omass = 15.999
Nmass = 14.007
Smass = 32.066

def usage():
    print __doc__

def add(x,y):
    return x+y

def sum(seq):
    return reduce(add, seq)

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 findCOM():
    #find center of mass
    Xcom = 0.0
    Ycom = 0.0
    Zcom = 0.0
    totalMass = 0.0

    for i in range(0,len(indices)):
        if (atypes[i] == "H"):
            myMass = Hmass
        elif (atypes[i] == "C"):
            myMass = Cmass
        elif (atypes[i] == "O"):
            myMass = Omass
        elif (atypes[i] == "N"):
            myMass = Nmass
        elif (atypes[i] == "S"):
            myMass = Smass
        else:
            print "unknown atom type!"
        
        Xcom = Xcom + myMass * positions[i][0]
        Ycom = Ycom + myMass * positions[i][1]
        Zcom = Zcom + myMass * positions[i][2]
        totalMass = totalMass + myMass

    Xcom = Xcom / totalMass
    Ycom = Ycom / totalMass
    Zcom = Zcom / totalMass

    COM = [Xcom, Ycom, Zcom]

    return COM

def findMoments():

    COM = findCOM()
    
    #find inertia tensor matrix elements

    I = numpy.zeros((3,3), numpy.float)
    
    for i in range(0,len(indices)):
        if (atypes[i] == "H"):
            myMass = Hmass
        elif (atypes[i] == "C"):
            myMass = Cmass
        elif (atypes[i] == "O"):
            myMass = Omass
        elif (atypes[i] == "N"):
            myMass = Nmass
        elif (atypes[i] == "S"):
            myMass = Smass
        else:
            print "unknown atom type!"

        dx = positions[i][0] - COM[0]
        dy = positions[i][1] - COM[1]
        dz = positions[i][2] - COM[2]

        I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
        I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
        I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )

        I[0,1] = I[0,1] - myMass * ( dx * dy )
        I[0,2] = I[0,2] - myMass * ( dx * dz )

        I[1,2] = I[1,2] - myMass * ( dy * dz )

        I[1,0] = I[0,1]
        I[2,0] = I[0,2]
        I[2,1] = I[1,2]

    print "Inertia Tensor:"
    print I
    print

    (evals, evects) = numpy.linalg.eig(I)
    print "evals:"
    print evals
    print 
    print "evects:"
    print evects
    print

    return (COM, evals, evects)

def writeFile(OutFileName):

    (COM, evals, evects) = findMoments()


    # we need to re-order the axes so that the smallest moment of inertia
    # (which corresponds to the long axis of the molecule) is along the z-axis
    # we'll just reverse the order of the three axes:
    
    axOrder = numpy.argsort(evals)    
    RotMat = numpy.zeros((3,3), numpy.float)
    RotMat[0] = evects[axOrder[2]]
    RotMat[1] = evects[axOrder[1]]
    RotMat[2] = evects[axOrder[0]]

    nAtoms = len(indices)
    
    print "writing output XYZ file"

    OutFile = open(OutFileName, 'w')        

    OutFile.write('%10d\n' % (nAtoms))
    OutFile.write('\n')

    for i in range(nAtoms):

        dx = positions[i][0] - COM[0]
        dy = positions[i][1] - COM[1]
        dz = positions[i][2] - COM[2]

        r = numpy.array([dx,dy,dz])
        rnew = numpy.dot(RotMat, r)
           
        OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
    OutFile.close()        

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="]) 
    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"): 
            XYZFileName = arg
            global _haveXYZFileName
            _haveXYZFileName = 1
        elif opt in ("-o", "--out"): 
            OutFileName = arg
            global _haveOutFileName
            _haveOutFileName = 1


    if (_haveXYZFileName != 1):
        usage() 
        print "No xyz file was specified"
        sys.exit()

    readFile(XYZFileName)

    if (_haveOutFileName == 1):
        writeFile(OutFileName)
    else:
        findMoments()
        
if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1442
Committed:	Mon May 10 17:28:26 2010 UTC (15 years, 5 months ago) by gezelter
File size:	6174 byte(s)
Log Message:	Adding property set to svn entries
#	User	Rev	Content
1	cpuglis	1189	#!/usr/bin/env python
2			"""principalAxisCalculator
3
4			Opens an XYZ file and computes the moments of inertia and principal axes
5	cli2	1360	for the structure in the XYZ file. Optionally rotate the structure so that
6			the long axis (that with the smallest eigenvalue) is pointing along the
7			z-axis.
8	cpuglis	1189
9			Usage: principalAxisCalculator
10
11			Options:
12			-h, --help show this help
13			-x, --xyz=... use specified XYZ (.xyz) file for the structure
14	cli2	1360	-o, --out=... rotate the structure so that the smallest eigenvalue
15			of the rotation matrix points along the z-axis.
16	cpuglis	1189
17			Example:
18	cli2	1360	principalAxisCalculator -x junk.xyz -o rot.xyz
19	cpuglis	1189
20			"""
21
22			__author__ = "Dan Gezelter (gezelter@nd.edu)"
23	gezelter	1442	__version__ = "$Revision$"
24			__date__ = "$Date$"
25	cpuglis	1189	__copyright__ = "Copyright (c) 2006 by the University of Notre Dame"
26	gezelter	1390	__license__ = "OpenMD"
27	cpuglis	1189
28			import sys
29			import getopt
30			import string
31			import math
32			import random
33			from sets import *
34	cli2	1360	import numpy
35	cpuglis	1189
36	cli2	1360
37	cpuglis	1189	_haveXYZFileName = 0
38	cli2	1360	_haveOutFileName = 0
39	cpuglis	1189
40			positions = []
41			indices = []
42			atypes = []
43	cli2	1360	Hmass = 1.0079
44			Cmass = 12.011
45			Omass = 15.999
46			Nmass = 14.007
47			Smass = 32.066
48	cpuglis	1189
49			def usage():
50			print __doc__
51
52			def add(x,y):
53			return x+y
54
55			def sum(seq):
56			return reduce(add, seq)
57
58			def readFile(XYZFileName):
59			print "reading XYZ file"
60
61			XYZFile = open(XYZFileName, 'r')
62			# Find number of atoms first
63			line = XYZFile.readline()
64			L = line.split()
65			nAtoms = int(L[0])
66			# skip comment line
67			line = XYZFile.readline()
68			for i in range(nAtoms):
69			line = XYZFile.readline()
70			L = line.split()
71			myIndex = i
72			indices.append(myIndex)
73			atomType = L[0]
74			atypes.append(atomType)
75			x = float(L[1])
76			y = float(L[2])
77			z = float(L[3])
78			positions.append([x, y, z])
79			XYZFile.close()
80
81
82	cli2	1360	def findCOM():
83	cpuglis	1189	#find center of mass
84			Xcom = 0.0
85			Ycom = 0.0
86			Zcom = 0.0
87			totalMass = 0.0
88
89	cli2	1360	for i in range(0,len(indices)):
90	cpuglis	1189	if (atypes[i] == "H"):
91			myMass = Hmass
92			elif (atypes[i] == "C"):
93			myMass = Cmass
94			elif (atypes[i] == "O"):
95			myMass = Omass
96			elif (atypes[i] == "N"):
97			myMass = Nmass
98			elif (atypes[i] == "S"):
99			myMass = Smass
100			else:
101			print "unknown atom type!"
102
103			Xcom = Xcom + myMass * positions[i][0]
104			Ycom = Ycom + myMass * positions[i][1]
105			Zcom = Zcom + myMass * positions[i][2]
106			totalMass = totalMass + myMass
107
108			Xcom = Xcom / totalMass
109			Ycom = Ycom / totalMass
110			Zcom = Zcom / totalMass
111
112			COM = [Xcom, Ycom, Zcom]
113
114	cli2	1360	return COM
115
116			def findMoments():
117
118			COM = findCOM()
119
120	cpuglis	1189	#find inertia tensor matrix elements
121
122	cli2	1360	I = numpy.zeros((3,3), numpy.float)
123	cpuglis	1189
124	cli2	1360	for i in range(0,len(indices)):
125	cpuglis	1189	if (atypes[i] == "H"):
126			myMass = Hmass
127			elif (atypes[i] == "C"):
128			myMass = Cmass
129			elif (atypes[i] == "O"):
130			myMass = Omass
131			elif (atypes[i] == "N"):
132			myMass = Nmass
133			elif (atypes[i] == "S"):
134			myMass = Smass
135			else:
136			print "unknown atom type!"
137
138	cli2	1360	dx = positions[i][0] - COM[0]
139			dy = positions[i][1] - COM[1]
140			dz = positions[i][2] - COM[2]
141	cpuglis	1189
142			I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
143			I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
144			I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )
145
146			I[0,1] = I[0,1] - myMass * ( dx * dy )
147			I[0,2] = I[0,2] - myMass * ( dx * dz )
148
149			I[1,2] = I[1,2] - myMass * ( dy * dz )
150
151			I[1,0] = I[0,1]
152			I[2,0] = I[0,2]
153			I[2,1] = I[1,2]
154
155			print "Inertia Tensor:"
156			print I
157			print
158
159	cli2	1360	(evals, evects) = numpy.linalg.eig(I)
160	cpuglis	1189	print "evals:"
161			print evals
162			print
163			print "evects:"
164			print evects
165			print
166
167	cli2	1360	return (COM, evals, evects)
168	cpuglis	1189
169	cli2	1360	def writeFile(OutFileName):
170	cpuglis	1189
171	cli2	1360	(COM, evals, evects) = findMoments()
172
173
174			# we need to re-order the axes so that the smallest moment of inertia
175			# (which corresponds to the long axis of the molecule) is along the z-axis
176			# we'll just reverse the order of the three axes:
177
178			axOrder = numpy.argsort(evals)
179			RotMat = numpy.zeros((3,3), numpy.float)
180			RotMat[0] = evects[axOrder[2]]
181			RotMat[1] = evects[axOrder[1]]
182			RotMat[2] = evects[axOrder[0]]
183
184			nAtoms = len(indices)
185
186			print "writing output XYZ file"
187
188			OutFile = open(OutFileName, 'w')
189
190			OutFile.write('%10d\n' % (nAtoms))
191			OutFile.write('\n')
192
193			for i in range(nAtoms):
194
195			dx = positions[i][0] - COM[0]
196			dy = positions[i][1] - COM[1]
197			dz = positions[i][2] - COM[2]
198
199			r = numpy.array([dx,dy,dz])
200			rnew = numpy.dot(RotMat, r)
201
202			OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
203			OutFile.close()
204
205	cpuglis	1189	def main(argv):
206			try:
207	cli2	1360	opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
208	cpuglis	1189	except getopt.GetoptError:
209			usage()
210			sys.exit(2)
211			for opt, arg in opts:
212			if opt in ("-h", "--help"):
213			usage()
214			sys.exit()
215			elif opt in ("-x", "--xyz"):
216			XYZFileName = arg
217			global _haveXYZFileName
218			_haveXYZFileName = 1
219	cli2	1360	elif opt in ("-o", "--out"):
220			OutFileName = arg
221			global _haveOutFileName
222			_haveOutFileName = 1
223	cpuglis	1189
224
225			if (_haveXYZFileName != 1):
226			usage()
227			print "No xyz file was specified"
228			sys.exit()
229
230			readFile(XYZFileName)
231
232	cli2	1360	if (_haveOutFileName == 1):
233			writeFile(OutFileName)
234			else:
235			findMoments()
236
237	cpuglis	1189	if __name__ == "__main__":
238			if len(sys.argv) == 1:
239			usage()
240			sys.exit()
241			main(sys.argv[1:])