[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]]

    q = [0.0, 0.0, 0.0, 0.0]
    myEuler = [0.0, 0.0, 0.0]

    # RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
    
    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        

    print "Quaternions:"
    print q

    theta = math.acos(RotMat[2][2])
    ctheta = RotMat[2][2]
    stheta = math.sqrt(1.0 - ctheta * ctheta)

    if (math.fabs(stheta) < 1e-6):
        psi = 0.0
        phi = math.atan2(-RotMat[1][0], RotMat[0][0]) 
    else:
        phi = math.atan2(RotMat[2][0], -RotMat[2][1])
        psi = math.atan2(RotMat[0][2], RotMat[1][2])
        
    if (phi < 0):
        phi = phi + 2.0 * math.pi;
        
    if (psi < 0):
        psi = psi + 2.0 * math.pi;
        
    myEuler[0] = phi * 180.0 / math.pi;
    myEuler[1] = theta * 180.0 / math.pi;
    myEuler[2] = psi * 180.0 / math.pi;

    print "Euler Angles:"
    print myEuler

    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:	1877
Committed:	Thu Jun 6 15:43:35 2013 UTC (12 years, 5 months ago) by gezelter
File size:	8477 byte(s)
Log Message:	New electrostatic method, starting to do some performance tuning.
#	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			# we need to re-order the axes so that the smallest moment of inertia
189			# (which corresponds to the long axis of the molecule) is along the z-axis
190			# we'll just reverse the order of the three axes:
191
192			axOrder = numpy.argsort(evals)
193			RotMat = numpy.zeros((3,3), numpy.float)
194			RotMat[0] = evects[axOrder[2]]
195			RotMat[1] = evects[axOrder[1]]
196			RotMat[2] = evects[axOrder[0]]
197
198	gezelter	1877	q = [0.0, 0.0, 0.0, 0.0]
199			myEuler = [0.0, 0.0, 0.0]
200
201			# RotMat to Quat code is out of OpenMD's SquareMatrix3.hpp code:
202
203			t = RotMat[0][0] + RotMat[1][1] + RotMat[2][2] + 1.0
204
205			if( t > 1e-6 ):
206			s = 0.5 / math.sqrt( t )
207			q[0] = 0.25 / s
208			q[1] = (RotMat[1][2] - RotMat[2][1]) * s
209			q[2] = (RotMat[2][0] - RotMat[0][2]) * s
210			q[3] = (RotMat[0][1] - RotMat[1][0]) * s
211			else:
212			ad1 = RotMat[0][0]
213			ad2 = RotMat[1][1]
214			ad3 = RotMat[2][2]
215
216			if( ad1 >= ad2 and ad1 >= ad3 ):
217			s = 0.5 / math.sqrt( 1.0 + RotMat[0][0] - RotMat[1][1] - RotMat[2][2] )
218			q[0] = (RotMat[1][2] - RotMat[2][1]) * s
219			q[1] = 0.25 / s
220			q[2] = (RotMat[0][1] + RotMat[1][0]) * s
221			q[3] = (RotMat[0][2] + RotMat[2][0]) * s
222			elif ( ad2 >= ad1 and ad2 >= ad3 ):
223			s = 0.5 / math.sqrt( 1.0 + RotMat[1][1] - RotMat[0][0] - RotMat[2][2] )
224			q[0] = (RotMat[2][0] - RotMat[0][2] ) * s
225			q[1] = (RotMat[0][1] + RotMat[1][0]) * s
226			q[2] = 0.25 / s
227			q[3] = (RotMat[1][2] + RotMat[2][1]) * s
228			else:
229			s = 0.5 / math.sqrt( 1.0 + RotMat[2][2] - RotMat[0][0] - RotMat[1][1] )
230			q[0] = (RotMat[0][1] - RotMat[1][0]) * s
231			q[1] = (RotMat[0][2] + RotMat[2][0]) * s
232			q[2] = (RotMat[1][2] + RotMat[2][1]) * s
233
234			print "Quaternions:"
235			print q
236
237			theta = math.acos(RotMat[2][2])
238			ctheta = RotMat[2][2]
239			stheta = math.sqrt(1.0 - ctheta * ctheta)
240
241			if (math.fabs(stheta) < 1e-6):
242			psi = 0.0
243			phi = math.atan2(-RotMat[1][0], RotMat[0][0])
244			else:
245			phi = math.atan2(RotMat[2][0], -RotMat[2][1])
246			psi = math.atan2(RotMat[0][2], RotMat[1][2])
247
248			if (phi < 0):
249			phi = phi + 2.0 * math.pi;
250
251			if (psi < 0):
252			psi = psi + 2.0 * math.pi;
253
254			myEuler[0] = phi * 180.0 / math.pi;
255			myEuler[1] = theta * 180.0 / math.pi;
256			myEuler[2] = psi * 180.0 / math.pi;
257
258			print "Euler Angles:"
259			print myEuler
260
261	gezelter	1798	nAtoms = len(indices)
262
263			print "writing output XYZ file"
264
265			OutFile = open(OutFileName, 'w')
266
267			OutFile.write('%10d\n' % (nAtoms))
268			OutFile.write('\n')
269
270			for i in range(nAtoms):
271
272			dx = positions[i][0] - COM[0]
273			dy = positions[i][1] - COM[1]
274			dz = positions[i][2] - COM[2]
275
276			r = numpy.array([dx,dy,dz])
277			rnew = numpy.dot(RotMat, r)
278
279			OutFile.write('%s\t%f\t%f\t%f\t%d\n' % (atypes[i], rnew[0], rnew[1], rnew[2], i))
280			OutFile.close()
281
282			def main(argv):
283			try:
284			opts, args = getopt.getopt(argv, "hx:o:", ["help", "xyz=", "out="])
285			except getopt.GetoptError:
286			usage()
287			sys.exit(2)
288			for opt, arg in opts:
289			if opt in ("-h", "--help"):
290			usage()
291			sys.exit()
292			elif opt in ("-x", "--xyz"):
293			XYZFileName = arg
294			global _haveXYZFileName
295			_haveXYZFileName = 1
296			elif opt in ("-o", "--out"):
297			OutFileName = arg
298			global _haveOutFileName
299			_haveOutFileName = 1
300
301
302			if (_haveXYZFileName != 1):
303			usage()
304			print "No xyz file was specified"
305			sys.exit()
306
307			readFile(XYZFileName)
308
309			if (_haveOutFileName == 1):
310			writeFile(OutFileName)
311			else:
312			findMoments()
313
314			if __name__ == "__main__":
315			if len(sys.argv) == 1:
316			usage()
317			sys.exit()
318			main(sys.argv[1:])