[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
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:	1798
Committed:	Thu Sep 13 14:10:11 2012 UTC (12 years, 7 months ago) by gezelter
File size:	6387 byte(s)
Log Message:	Merged trunk changes into the development branch
#	User	Rev	Content
1	gezelter	1798	#!@PYTHON_EXECUTABLE@
2			"""principalAxisCalculator
3
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	kstocke1	1701	Aumass = 196.466569
49	gezelter	1528	Au1mass = 196.466569
50			Au2mass = 0.5
51	gezelter	1798
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			def findCOM():
86			#find center of mass
87			Xcom = 0.0
88			Ycom = 0.0
89			Zcom = 0.0
90			totalMass = 0.0
91
92			for i in range(0,len(indices)):
93			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			elif (atypes[i] == "Au1"):
104			myMass = Au1mass
105			elif (atypes[i] == "Au2"):
106			myMass = Au2mass
107			elif (atypes[i] == "Au"):
108			myMass = Aumass
109			else:
110			print "unknown atom type! %s" % (atypes[i])
111
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			return COM
124
125			def findMoments():
126
127			COM = findCOM()
128
129			#find inertia tensor matrix elements
130
131			I = numpy.zeros((3,3), numpy.float)
132
133			for i in range(0,len(indices)):
134			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			elif (atypes[i] == "Au1"):
145			myMass = Au1mass
146			elif (atypes[i] == "Au2"):
147			myMass = Au2mass
148			elif (atypes[i] == "Au"):
149			myMass = Aumass
150			else:
151			print "unknown atom type!"
152
153			dx = positions[i][0] - COM[0]
154			dy = positions[i][1] - COM[1]
155			dz = positions[i][2] - COM[2]
156
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			(evals, evects) = numpy.linalg.eig(I)
175			print "evals:"
176			print evals
177			print
178			print "evects:"
179			print evects
180			print
181
182			return (COM, evals, evects)
183
184			def writeFile(OutFileName):
185
186			(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			def main(argv):
221			try:
222			opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
223			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			elif opt in ("-o", "--out"):
235			OutFileName = arg
236			global _haveOutFileName
237			_haveOutFileName = 1
238
239
240			if (_haveXYZFileName != 1):
241			usage()
242			print "No xyz file was specified"
243			sys.exit()
244
245			readFile(XYZFileName)
246
247			if (_haveOutFileName == 1):
248			writeFile(OutFileName)
249			else:
250			findMoments()
251
252			if __name__ == "__main__":
253			if len(sys.argv) == 1:
254			usage()
255			sys.exit()
256			main(sys.argv[1:])