[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
Aumass = 196.466569
Au1mass = 196.466569
Au2mass =   0.5

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
        elif (atypes[i] == "Au1"):
            myMass = Au1mass
        elif (atypes[i] == "Au2"):
            myMass = Au2mass
        elif (atypes[i] == "Au"):
            myMass = Aumass
        else:
            print "unknown atom type! %s" % (atypes[i])
        
        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
        elif (atypes[i] == "Au1"):
            myMass = Au1mass
        elif (atypes[i] == "Au2"):
            myMass = Au2mass
        elif (atypes[i] == "Au"):
            myMass = Aumass
        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:	1782
Committed:	Wed Aug 22 02:28:28 2012 UTC (12 years, 10 months ago) by gezelter
File size:	6640 byte(s)
Log Message:	MERGE OpenMD development branch 1465:1781 into trunk
#	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	gezelter	1782	Aumass = 196.466569
49			Au1mass = 196.466569
50			Au2mass = 0.5
51	cpuglis	1189
52			def usage():
53			print __doc__
54
55			def add(x,y):
56			return x+y
57
58			def sum(seq):
59			return reduce(add, seq)
60
61			def readFile(XYZFileName):
62			print "reading XYZ file"
63
64			XYZFile = open(XYZFileName, 'r')
65			# Find number of atoms first
66			line = XYZFile.readline()
67			L = line.split()
68			nAtoms = int(L[0])
69			# skip comment line
70			line = XYZFile.readline()
71			for i in range(nAtoms):
72			line = XYZFile.readline()
73			L = line.split()
74			myIndex = i
75			indices.append(myIndex)
76			atomType = L[0]
77			atypes.append(atomType)
78			x = float(L[1])
79			y = float(L[2])
80			z = float(L[3])
81			positions.append([x, y, z])
82			XYZFile.close()
83
84
85	cli2	1360	def findCOM():
86	cpuglis	1189	#find center of mass
87			Xcom = 0.0
88			Ycom = 0.0
89			Zcom = 0.0
90			totalMass = 0.0
91
92	cli2	1360	for i in range(0,len(indices)):
93	cpuglis	1189	if (atypes[i] == "H"):
94			myMass = Hmass
95			elif (atypes[i] == "C"):
96			myMass = Cmass
97			elif (atypes[i] == "O"):
98			myMass = Omass
99			elif (atypes[i] == "N"):
100			myMass = Nmass
101			elif (atypes[i] == "S"):
102			myMass = Smass
103	gezelter	1782	elif (atypes[i] == "Au1"):
104			myMass = Au1mass
105			elif (atypes[i] == "Au2"):
106			myMass = Au2mass
107			elif (atypes[i] == "Au"):
108			myMass = Aumass
109	cpuglis	1189	else:
110	gezelter	1782	print "unknown atom type! %s" % (atypes[i])
111	cpuglis	1189
112			Xcom = Xcom + myMass * positions[i][0]
113			Ycom = Ycom + myMass * positions[i][1]
114			Zcom = Zcom + myMass * positions[i][2]
115			totalMass = totalMass + myMass
116
117			Xcom = Xcom / totalMass
118			Ycom = Ycom / totalMass
119			Zcom = Zcom / totalMass
120
121			COM = [Xcom, Ycom, Zcom]
122
123	cli2	1360	return COM
124
125			def findMoments():
126
127			COM = findCOM()
128
129	cpuglis	1189	#find inertia tensor matrix elements
130
131	cli2	1360	I = numpy.zeros((3,3), numpy.float)
132	cpuglis	1189
133	cli2	1360	for i in range(0,len(indices)):
134	cpuglis	1189	if (atypes[i] == "H"):
135			myMass = Hmass
136			elif (atypes[i] == "C"):
137			myMass = Cmass
138			elif (atypes[i] == "O"):
139			myMass = Omass
140			elif (atypes[i] == "N"):
141			myMass = Nmass
142			elif (atypes[i] == "S"):
143			myMass = Smass
144	gezelter	1782	elif (atypes[i] == "Au1"):
145			myMass = Au1mass
146			elif (atypes[i] == "Au2"):
147			myMass = Au2mass
148			elif (atypes[i] == "Au"):
149			myMass = Aumass
150	cpuglis	1189	else:
151			print "unknown atom type!"
152
153	cli2	1360	dx = positions[i][0] - COM[0]
154			dy = positions[i][1] - COM[1]
155			dz = positions[i][2] - COM[2]
156	cpuglis	1189
157			I[0,0] = I[0,0] + myMass * ( dy * dy + dz * dz )
158			I[1,1] = I[1,1] + myMass * ( dx * dx + dz * dz )
159			I[2,2] = I[2,2] + myMass * ( dx * dx + dy * dy )
160
161			I[0,1] = I[0,1] - myMass * ( dx * dy )
162			I[0,2] = I[0,2] - myMass * ( dx * dz )
163
164			I[1,2] = I[1,2] - myMass * ( dy * dz )
165
166			I[1,0] = I[0,1]
167			I[2,0] = I[0,2]
168			I[2,1] = I[1,2]
169
170			print "Inertia Tensor:"
171			print I
172			print
173
174	cli2	1360	(evals, evects) = numpy.linalg.eig(I)
175	cpuglis	1189	print "evals:"
176			print evals
177			print
178			print "evects:"
179			print evects
180			print
181
182	cli2	1360	return (COM, evals, evects)
183	cpuglis	1189
184	cli2	1360	def writeFile(OutFileName):
185	cpuglis	1189
186	cli2	1360	(COM, evals, evects) = findMoments()
187
188
189			# we need to re-order the axes so that the smallest moment of inertia
190			# (which corresponds to the long axis of the molecule) is along the z-axis
191			# we'll just reverse the order of the three axes:
192
193			axOrder = numpy.argsort(evals)
194			RotMat = numpy.zeros((3,3), numpy.float)
195			RotMat[0] = evects[axOrder[2]]
196			RotMat[1] = evects[axOrder[1]]
197			RotMat[2] = evects[axOrder[0]]
198
199			nAtoms = len(indices)
200
201			print "writing output XYZ file"
202
203			OutFile = open(OutFileName, 'w')
204
205			OutFile.write('%10d\n' % (nAtoms))
206			OutFile.write('\n')
207
208			for i in range(nAtoms):
209
210			dx = positions[i][0] - COM[0]
211			dy = positions[i][1] - COM[1]
212			dz = positions[i][2] - COM[2]
213
214			r = numpy.array([dx,dy,dz])
215			rnew = numpy.dot(RotMat, r)
216
217			OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
218			OutFile.close()
219
220	cpuglis	1189	def main(argv):
221			try:
222	cli2	1360	opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
223	cpuglis	1189	except getopt.GetoptError:
224			usage()
225			sys.exit(2)
226			for opt, arg in opts:
227			if opt in ("-h", "--help"):
228			usage()
229			sys.exit()
230			elif opt in ("-x", "--xyz"):
231			XYZFileName = arg
232			global _haveXYZFileName
233			_haveXYZFileName = 1
234	cli2	1360	elif opt in ("-o", "--out"):
235			OutFileName = arg
236			global _haveOutFileName
237			_haveOutFileName = 1
238	cpuglis	1189
239
240			if (_haveXYZFileName != 1):
241			usage()
242			print "No xyz file was specified"
243			sys.exit()
244
245			readFile(XYZFileName)
246
247	cli2	1360	if (_haveOutFileName == 1):
248			writeFile(OutFileName)
249			else:
250			findMoments()
251
252	cpuglis	1189	if __name__ == "__main__":
253			if len(sys.argv) == 1:
254			usage()
255			sys.exit()
256			main(sys.argv[1:])