[ViewVC] Contents of: OpenMD/trunk/src/applications/hydrodynamics/diffExplainer

#!@PYTHON_EXECUTABLE@
"""Computes predicted diffusion constants and rotational relaxation
times from a diff file.  Explains the values in the diff file in
terms of properties that can be calculated from a molecular dynamics
simulation.

Usage: diffExplainer

Options:
  -h, --help              show this help
  -d, --diff-file=...     use specified diff file
  -t, --temperature=...   temperature in Kelvin


Example:
   diffExplainer -d ellipsoid.diff -t 300

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision$"
__date__ = "$Date$"

__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math

_haveDiffFileName = 0
_haveTemperature = 0

def usage():
    print __doc__

def readDiffFile(diffFileName):
    diffFile = open(diffFileName, 'r')

    global DTTxx
    global DTTyy
    global DTTzz

    global XiRRxx
    global XiRRyy
    global XiRRzz

    global DRRxx
    global DRRyy
    global DRRzz
     
    line = diffFile.readline()
    while 1:
        L = line.split()
        XiRRxx = float(L[31])
        XiRRyy = float(L[35])
        XiRRzz = float(L[39])
        DTTxx = float(L[115])
        DTTyy = float(L[119])
        DTTzz = float(L[123])
        DRRxx = float(L[142])
        DRRyy = float(L[146])
        DRRzz = float(L[150])              
        line = diffFile.readline()
        if not line: break
        
    diffFile.close()
    
def computeProperties(temperature):
    
    kb = 1.9872156E-3    
    kt = kb * temperature

    print
    print "Translational Diffusion Constant (angstroms^2 fs^-1): %.3e" % ((DTTxx + DTTyy + DTTzz)/3.0)

    Delta = math.sqrt(math.pow((DRRxx-DRRyy),2) +(DRRzz-DRRxx)*(DRRzz-DRRyy))
    a = math.sqrt(3.0)*(DRRxx-DRRyy) 
    b = (2.0*DRRzz - DRRxx - DRRyy + 2.0*Delta)
    N = 2.0*math.sqrt(Delta)*math.sqrt(b)
    Di = (DRRxx + DRRyy + DRRzz)/3.0
    f1 = 6.0*Di + 2.0*Delta
    f2 = 6.0*Di - 2.0*Delta
    f3 = 3.0*(DRRxx + Di)
    f4 = 3.0*(DRRyy + Di)
    f5 = 3.0*(DRRzz + Di)

    f0 = (f1 + f2 + f3 + f4 + f5)/5.0

    print
    print "l=2 Orientational correlation functions:"
    print
    print "In general only the C^2(0,0) correlation function relates to the lcorr"
    print "computed via the DynamicProps correlation function routine."
    print 
    print "To get any of the specific correlation functions, multiply each amplitude"
    print "by exp(-t/t_i) where t_i is the relevant decay time printed in the first row:"
    print
    print "decay times (fs):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (1.0/f1, 1.0/f2, 1.0/f3, 1.0/f4, 1.0/f5 )
    print
    print "       C^2(0, 0):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (pow(a/N,2), pow(b/N,2), 0, 0, 0 )
    print "       C^2(1, 1):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (0,0,0.5,0.5,0)
    print "       C^2(1,-1):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (0,0,-0.5,0.5,0)
    print "       C^2(2, 2):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (0.5*pow(b/N,2), 0.5*pow(a/N,2), 0, 0, 0.5)
    print "       C^2(2,-2):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (0.5*pow(b/N,2), 0.5*pow(a/N,2), 0, 0,-0.5)
    print "       C^2(2, 0):   %.3e   %.3e   %.3e   %.3e   %.3e"  % (math.sqrt(2.0)*a*b/pow(N,2),-math.sqrt(2.0)*a*b/pow(N,2),0,0,0)
    print
    print 
    print "average (or characteristic) relaxation time:\t%.3e" % (1.0/f0)

def main(argv):                         
    try:                                
        opts, args = getopt.getopt(argv, "hd:t:", ["help", "diff-file=", "temperature="]) 
    except getopt.GetoptError:           
        usage()                          
        sys.exit(2)                     
    for opt, arg in opts:                
        if opt in ("-h", "--help"):      
            usage()                     
            sys.exit()                  
        elif opt in ("-d", "--diff-file"): 
            diffFileName = arg
            global _haveDiffFileName
            _haveDiffFileName = 1
        elif opt in ("-t", "--temperature"): 
            temperature = float(arg)
            global _haveTemperature
            _haveTemperature = 1
    if (_haveDiffFileName != 1):
        usage() 
        print "No diff file was specified"
        sys.exit()
    if (_haveTemperature != 1):
        usage() 
        print "No temperature was specified"
        sys.exit()
        
    readDiffFile(diffFileName);
    computeProperties(temperature);

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1795
Committed:	Fri Sep 7 18:13:55 2012 UTC (13 years, 1 month ago) by gezelter
File size:	4518 byte(s)
Log Message:	Windows fixes
#	Content
1	#!@PYTHON_EXECUTABLE@
2	"""Computes predicted diffusion constants and rotational relaxation
3	times from a diff file. Explains the values in the diff file in
4	terms of properties that can be calculated from a molecular dynamics
5	simulation.
6
7	Usage: diffExplainer
8
9	Options:
10	-h, --help show this help
11	-d, --diff-file=... use specified diff file
12	-t, --temperature=... temperature in Kelvin
13
14
15	Example:
16	diffExplainer -d ellipsoid.diff -t 300
17
18	"""
19
20	__author__ = "Dan Gezelter (gezelter@nd.edu)"
21	__version__ = "$Revision$"
22	__date__ = "$Date$"
23
24	__copyright__ = "Copyright (c) 2008 by the University of Notre Dame"
25	__license__ = "OpenMD"
26
27	import sys
28	import getopt
29	import string
30	import math
31
32	_haveDiffFileName = 0
33	_haveTemperature = 0
34
35	def usage():
36	print __doc__
37
38	def readDiffFile(diffFileName):
39	diffFile = open(diffFileName, 'r')
40
41	global DTTxx
42	global DTTyy
43	global DTTzz
44
45	global XiRRxx
46	global XiRRyy
47	global XiRRzz
48
49	global DRRxx
50	global DRRyy
51	global DRRzz
52
53	line = diffFile.readline()
54	while 1:
55	L = line.split()
56	XiRRxx = float(L[31])
57	XiRRyy = float(L[35])
58	XiRRzz = float(L[39])
59	DTTxx = float(L[115])
60	DTTyy = float(L[119])
61	DTTzz = float(L[123])
62	DRRxx = float(L[142])
63	DRRyy = float(L[146])
64	DRRzz = float(L[150])
65	line = diffFile.readline()
66	if not line: break
67
68	diffFile.close()
69
70	def computeProperties(temperature):
71
72	kb = 1.9872156E-3
73	kt = kb * temperature
74
75	print
76	print "Translational Diffusion Constant (angstroms^2 fs^-1): %.3e" % ((DTTxx + DTTyy + DTTzz)/3.0)
77
78	Delta = math.sqrt(math.pow((DRRxx-DRRyy),2) +(DRRzz-DRRxx)*(DRRzz-DRRyy))
79	a = math.sqrt(3.0)*(DRRxx-DRRyy)
80	b = (2.0DRRzz - DRRxx - DRRyy + 2.0Delta)
81	N = 2.0math.sqrt(Delta)math.sqrt(b)
82	Di = (DRRxx + DRRyy + DRRzz)/3.0
83	f1 = 6.0Di + 2.0Delta
84	f2 = 6.0Di - 2.0Delta
85	f3 = 3.0*(DRRxx + Di)
86	f4 = 3.0*(DRRyy + Di)
87	f5 = 3.0*(DRRzz + Di)
88
89	f0 = (f1 + f2 + f3 + f4 + f5)/5.0
90
91	print
92	print "l=2 Orientational correlation functions:"
93	print
94	print "In general only the C^2(0,0) correlation function relates to the lcorr"
95	print "computed via the DynamicProps correlation function routine."
96	print
97	print "To get any of the specific correlation functions, multiply each amplitude"
98	print "by exp(-t/t_i) where t_i is the relevant decay time printed in the first row:"
99	print
100	print "decay times (fs): %.3e %.3e %.3e %.3e %.3e" % (1.0/f1, 1.0/f2, 1.0/f3, 1.0/f4, 1.0/f5 )
101	print
102	print " C^2(0, 0): %.3e %.3e %.3e %.3e %.3e" % (pow(a/N,2), pow(b/N,2), 0, 0, 0 )
103	print " C^2(1, 1): %.3e %.3e %.3e %.3e %.3e" % (0,0,0.5,0.5,0)
104	print " C^2(1,-1): %.3e %.3e %.3e %.3e %.3e" % (0,0,-0.5,0.5,0)
105	print " C^2(2, 2): %.3e %.3e %.3e %.3e %.3e" % (0.5pow(b/N,2), 0.5pow(a/N,2), 0, 0, 0.5)
106	print " C^2(2,-2): %.3e %.3e %.3e %.3e %.3e" % (0.5pow(b/N,2), 0.5pow(a/N,2), 0, 0,-0.5)
107	print " C^2(2, 0): %.3e %.3e %.3e %.3e %.3e" % (math.sqrt(2.0)ab/pow(N,2),-math.sqrt(2.0)ab/pow(N,2),0,0,0)
108	print
109	print
110	print "average (or characteristic) relaxation time:\t%.3e" % (1.0/f0)
111
112	def main(argv):
113	try:
114	opts, args = getopt.getopt(argv, "hd:t:", ["help", "diff-file=", "temperature="])
115	except getopt.GetoptError:
116	usage()
117	sys.exit(2)
118	for opt, arg in opts:
119	if opt in ("-h", "--help"):
120	usage()
121	sys.exit()
122	elif opt in ("-d", "--diff-file"):
123	diffFileName = arg
124	global _haveDiffFileName
125	_haveDiffFileName = 1
126	elif opt in ("-t", "--temperature"):
127	temperature = float(arg)
128	global _haveTemperature
129	_haveTemperature = 1
130	if (_haveDiffFileName != 1):
131	usage()
132	print "No diff file was specified"
133	sys.exit()
134	if (_haveTemperature != 1):
135	usage()
136	print "No temperature was specified"
137	sys.exit()
138
139	readDiffFile(diffFileName);
140	computeProperties(temperature);
141
142	if __name__ == "__main__":
143	if len(sys.argv) == 1:
144	usage()
145	sys.exit()
146	main(sys.argv[1:])