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

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