OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
solLiqFricCalc
1#!@Python3_EXECUTABLE@
2
3__author__ = "Patrick Louden (plouden@nd.edu)"
4__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
5__license__ = "OpenMD"
6
7import sys
8import os
9import datetime
10import argparse
11import textwrap
12import numpy as np
13import math
14from scipy.optimize import curve_fit
15from argparse import RawDescriptionHelpFormatter
16from scipy import stats
17
18def usage():
19 print(__doc__)
20
21# Fit function we are optimizing data to (using tanh)
22def funcTanh(x, Vs, Z1, Z2, ml, w, Vl):
23 ''' from Mathematica CForm
24 0.5*(Vl + ml*(x - Z1))*(1 - Tanh((x - Z1)/w)) +
25 0.5*Vs*(Tanh((x - Z1)/w) - Tanh((x - Z2)/w)) +
26 0.5*(Vl - ml*(x - Z2))*(1 + Tanh((x - Z2)/w))
27 '''
28
29 return ( 0.5*(Vl+ml*(x-Z1))*(1.0-np.tanh((x-Z1)/w))+0.5*Vs*(np.tanh((x-Z1)/w)-np.tanh((x-Z2)/w))+0.5*(Vl-ml*(x-Z2))*(1.0+np.tanh((x-Z2)/w)) )
30
31# Fit function we are optimizing data to (using piecewise smooth parabolic)
32def funcParab(x, Vs, Z1, Z2, ml, w, Vl):
33
34 y = np.empty_like (x)
35
36 k = 2.0 * (Vs - Vl - ml * (Z1-w)) / (w*w)
37
38 i = 0
39 for xi in x:
40 if (xi < (Z1-w)):
41 yi = Vl + ml * xi
42 elif ((xi > Z1-w) and (xi < Z1) ):
43 yi = Vs - k * (xi - Z1)*(xi - Z1) / 2.0
44 elif ((xi > Z1) and (xi < Z2)):
45 yi = Vs
46 elif ((xi > Z2) and (xi < Z2 + w)):
47 yi = Vs - k * (xi - Z2)*(xi - Z2) / 2.0
48 else:
49 yi = Vs - k * w * w / 2.0 - ml * (xi - Z2 -w)
50 y[i] = yi
51 i += 1
52
53 return y
54
55def funcParabOnePoint(x, Vs, Z1, Z2, ml, w, Vl):
56
57 k = 2.0 * (Vs - Vl - ml * (Z1-w)) / (w*w)
58 if (x < (Z1-w)):
59 y = Vl + ml * x
60 elif ((x > Z1-w) and (x < Z1) ):
61 y = Vs - k * (x - Z1)*(x - Z1) / 2.0
62 elif ((x > Z1) and (x < Z2)):
63 y = Vs
64 elif ((x > Z2) and (x < Z2 + w)):
65 y = Vs - k * (x - Z2)*(x - Z2) / 2.0
66 else:
67 y = Vs - k * w * w / 2.0 - ml * (x - Z2 -w)
68 return y
69
70# This function reads in the .rnemd file and extracts the actual
71# momentum flux, as well as the zPositions and the average velocities
72# therein.
73def rnemdExtractor(rnemdFileName):
74 if os.path.exists(rnemdFileName):
75 rnemdFile = open(rnemdFileName, 'r')
76 else:
77 #print 'Error: cannot open ' + rnemdFileName
78 sys.exit(2)
79
80 rnemdZPos = np.array([])
81 rnemdV = np.array([]).reshape(0, 3)
82
83 while True:
84 line = rnemdFile.readline()
85 if not line: break
86 if "Actual flux:" in line or "actual flux:" in line:
87 rnemdFile.readline()
88 line = rnemdFile.readline()
89 rnemdJzP = np.array([[float(line.split()[4][:-1])], [float(line.split()[5][:-1])], [float(line.split()[6])]])
90 elif "#" not in line:
91 if len(line.split()) > 6:
92 print("Error in the .rnemd file, found more than 6 expected columns.")
93 else:
94 rnemdZPos = np.append(rnemdZPos, np.array([float(line.split()[0])]))
95 rnemdV = np.vstack( [rnemdV, [ float(line.split()[2]), float(line.split()[3]), float(line.split()[4]) ] ] )
96
97 # Now let's determine what kind of shearing simulation was performed
98 # ie, JzPx, JzPy, or JzPxy
99 if (rnemdJzP[2] != 0.0):
100 print("Error: JzPz Actual Exchange total was found to be >0.0 in (.rnemd) file.")
101 sys.exit(2)
102 if (rnemdJzP[0] == 0.0 and rnemdJzP[1] == 0.0):
103 print("Error: both JzPx and JzPy Actual Exchange totals were found to be 0.0 in (.rnemd) file. ")
104 sys.exit(2)
105 elif (rnemdJzP[0] != 0.0 and rnemdJzP[1] != 0.0):
106 print("Error: both JzPx and JzPy Actual Exchange totals were found to be >0.0 in (.rnemd) file.\n This script only calculates friction coefficients for single dimensional momentum fluxes.")
107 sys.exit(2)
108
109 # After sanity checking, we now set a variable jzType to store if the
110 # imposed momentum flux is in the x- (0) or y-direction (1).
111 # jzType will be used when we iterate through rnemdV
112 # jzType = 0 -> flux was in x-dimension
113 # jzType = 1 -> flux was in y-dimension
114
115 elif (rnemdJzP[0] == 0.0 or rnemdJzP[1] == 0.0):
116 if np.abs(rnemdJzP[0]) > np.abs(rnemdJzP [1]):
117 jzType = 0
118 else:
119 jzType = 1
120
121 # print "Read in " + str(rnemdFileName)
122 return (rnemdJzP, jzType, rnemdZPos, rnemdV)
123
124
125def velFitter(rnemdJzP, jzType, rnemdZPos, rnemdV, Vs, Z1, Z2, m, w, Vl, numDel, outputFileName, useZLocations):
126 rnemdVx = np.array([]).reshape(0, 1)
127 rnemdVy = np.array([]).reshape(0, 1)
128
129 for i in range(0, len(rnemdV)):
130 rnemdVx = np.append(rnemdVx, rnemdV[i][0])
131 rnemdVy = np.append(rnemdVy, rnemdV[i][1])
132
133 # Remove the first two and last two values of rnemdVx and rnemdVy to smooth fits
134 for i in range(0, numDel):
135 rnemdVx = np.delete(rnemdVx, 0)
136 rnemdVy = np.delete(rnemdVy, 0)
137 rnemdZPos = np.delete(rnemdZPos, 0)
138 for i in range(0, numDel+1):
139 rnemdVx = np.delete(rnemdVx, len(rnemdVx)-1)
140 rnemdVy = np.delete(rnemdVy, len(rnemdVy)-1)
141 rnemdZPos = np.delete(rnemdZPos, len(rnemdZPos)-1)
142
143 # Tanh fitting function...
144 # curve_fit(function to fit to, x-data, y-data, [initial guess of param 1, initial guess of param 2,...]
145 if (jzType == 0):
146 popt, pcov = curve_fit(funcTanh, rnemdZPos, rnemdVx, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
147 elif (jzType == 1):
148 popt, pcov = curve_fit(funcTanh, rnemdZPos, rnemdVy, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
149
150 # Parabola fitting function...
151 # curve_fit(function to fit to, x-data, y-data, [initial guess of param 1, initial guess of param 2,...]
152
153 param_bounds = ([-np.inf, 0.0, 0.0, -np.inf, 0.0, -np.inf], [np.inf, np.inf, np.inf, np.inf, np.inf, np.inf])
154 if (jzType == 0):
155 poptParab, pcovParab = curve_fit(funcParab, rnemdZPos, rnemdVx, [Vs, Z1, Z2, m, w, Vl], bounds=param_bounds, ftol=0.01)
156 elif (jzType == 1):
157 poptParab, pcovParab = curve_fit(funcParab, rnemdZPos, rnemdVy, [Vs, Z1, Z2, m, w, Vl], bounds=param_bounds, ftol=0.01)
158
159 # write out the raw data and the fit for visual inspection of quality of fit
160 # print "Writing out fit file"
161 fitFileName = outputFileName + "FIT"
162 fitFile = open(fitFileName, "w")
163 fitFile.write("# zPosition <V_{i}> fit value" + "\n")
164
165 #y = np.empty_like(rnemdZPos)
166 #y = funcParab(rnemdZPos, poptParab[0], poptParab[1], poptParab[2], poptParab[3], poptParab[4], poptParab[5])
167
168 for i in range(0, len(rnemdZPos)):
169 y = funcParabOnePoint(rnemdZPos[i], poptParab[0], poptParab[1], poptParab[2], poptParab[3], poptParab[4], poptParab[5])
170 if (jzType == 0):
171 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVx[i]) + "\t" + str(y) + "\n")
172 elif (jzType == 1):
173 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVy[i]) + "\t" + str(y) + "\n")
174
175 # print "popt = ", popt
176 return (popt, pcov, poptParab, pcovParab)
177
178
179def calcShearRate(jzType, rnemdZPos, rnemdV):
180 # jzType = 0 -> flux was in x-dimension
181 # jzType = 1 -> flux was in y-dimension
182
183 velLiq = (rnemdV[0][jzType] + rnemdV[len(rnemdV)-1][jzType])/2.0
184
185 if (len(rnemdV) % 2) == 0 :
186 velSol = ( rnemdV[int(len(rnemdV) / 2)-1][jzType] + rnemdV[int(len(rnemdV)/ 2)][jzType]) / 2.0
187 elif (len(rnemdV) % 2) != 0:
188 velSol = ( rnemdV[int((len(rnemdV) + 1)/2)][jzType] + rnemdV[int((len(rnemdV)+1)/2)+1][jzType]) / 2.0
189
190 shearRate = velSol - velLiq
191
192 return (shearRate)
193
194
195def lambdaCalc(rnemdJzP, jzType, popt, sigma, convFactor):
196 # print "Performing lambda calculation"
197 # popt = Vs, Z1, Z2, ml, w, Vl
198
199 zLiq = (popt[0] - popt[5]) / popt[3] + popt[1]
200 if (popt[4]*convFactor < sigma):
201 zSol = popt[1] + 0.5*sigma
202 else:
203 zSol = popt[1] + 0.5*(popt[4]*convFactor)
204 delta = zLiq - zSol
205
206 if (jzType == 0):
207 lambda1 = rnemdJzP[0] / (delta * popt[3])
208 elif (jzType == 1):
209 lambda1 = rnemdJzP[1] / (delta * popt[3])
210
211 return (zLiq, zSol, delta, lambda1)
212
213
214def kappaCalc(rnemdJzP, jzType, z1, z2, w, popt, sigma, convFactor, useWidth, useZLocations):
215 # print "Performing kappa calculation"
216 # print "useZLocations = " , useZLocations
217
218 if (useZLocations):
219 print("Using supplied locations of the interface for kappa")
220 if (useWidth):
221 zLiq1 = z1 - 0.5*w
222 zSol1 = z1 + 0.5*w
223 zLiq2 = z2 + 0.5*w
224 zSol2 = z2 - 0.5*w
225 elif (not useWidth):
226 zLiq1 = z1 - 1.5*sigma
227 zSol1 = z1 + 1.5*sigma
228 zLiq2 = z2 + 1.5*sigma
229 zSol2 = z2 - 1.5*sigma
230 elif (not useZLocations):
231 print("Using the fit values for the locations of the interface for kappa")
232 if (useWidth):
233 zLiq1 = popt[1] - 0.5*w
234 zSol1 = popt[1] + 0.5*w
235 zLiq2 = popt[2] + 0.5*w
236 zSol2 = popt[2] - 0.5*w
237 elif (not useWidth):
238 zLiq1 = popt[1] - 1.5*sigma
239 zSol1 = popt[1] + 1.5*sigma
240 zLiq2 = popt[2] + 1.5*sigma
241 zSol2 = popt[2] - 1.5*sigma
242
243 vLiq1 = funcParabOnePoint(zLiq1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
244 vSol1 = funcParabOnePoint(zSol1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
245 deltaV1 = vSol1 - vLiq1
246
247 vLiq2 = funcParabOnePoint(zLiq2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
248 vSol2 = funcParabOnePoint(zSol2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
249 deltaV2 = vSol2 - vLiq2
250
251 if (jzType == 0):
252 k1 = rnemdJzP[0] / deltaV1
253 k2 = rnemdJzP[0] / deltaV2
254 elif (jzType == 1):
255 k1 = rnemdJzP[1] / deltaV1
256 k2 = rnemdJzP[1] / deltaV2
257
258 return (zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, k1, k2)
259
260
261def writeOutputFile(outputFileName, rnemdFileName, popt):
262 outFile = open(outputFileName, "w")
263
264 outFile.write("##################################################### \n")
265 outFile.write("## This output file was generated by solLiqFricCalc on " + str(datetime.datetime.now()) + " ## \n#\n")
266 outFile.write("# The velocity profile found in " + str(rnemdFileName) + " was fit by \n")
267 outFile.write("# \ty = Vl - ml * x for 0 <= x < (Z1 - w)\n")
268 outFile.write("# \ty = Vs - 0.5 * k * (x - Z1)^2 for (Z1 - w) <= x < Z1\n")
269 outFile.write("# \ty = Vs for Z1 <= x < Z2\n")
270 outFile.write("# \ty = Vs - 0.5 * k * (x - Z2)^2 for Z2 <= x < ( Z2 + w)\n")
271 outFile.write("# \ty = Vs - 0.5*k*w^2 - ml*(x-(Z2 + w)) for (Z2 + w) <= x\n")
272 outFile.write("# where (Vs, Vl, Z1, Z2, ml, w) are fit parameters, and here y = V(i) and x = z \n")
273 outFile.write("# Vs = velocity of the solid\n# Vl = velocity of the liquid at the interface \n")
274 outFile.write("# Z1 = z-position of the lower interface \n# Z2 = z-position of the upper interface \n")
275 outFile.write("# ml = slope of the velocity profile in the lower liquid region of the box \n")
276
277 outFile.write("# w = width of the interface \n#\n")
278 outFile.write("# Obtained optimized parameters \n")
279 outFile.write("# \t Vs = " + str(popt[0]) + "\n")
280 outFile.write("# \t Vl = " + str(popt[5]) + "\n")
281 outFile.write("# \t Z1 = " + str(popt[1]) + "\n")
282 outFile.write("# \t Z2 = " + str(popt[2]) + "\n")
283 outFile.write("# \t ml = " + str(popt[3]) + "\n")
284 outFile.write("# \t w = " + str(popt[4]) + "\n#\n")
285
286 if (float(popt[0]) < float(popt[5])):
287 print("Bad fit detected, Vs < Vl")
288 outFile.write("# WARNING Vs < Vl \n#\n")
289
290
291def writeOutputLambda(outputFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate):
292 outFile = open(outputFileName, "a")
293 outFile.write("# The interfacial friction coefficient, lambda, was calculated by \n")
294 outFile.write("# \tlambda = Jz(p) / (dV/dz)_{liquid} / delta \n")
295 outFile.write("# where Jz(p) is the imposed momentum flux, (dV/dz)_{liquid} is the slope of the velocity profile in the liquid \n")
296 outFile.write("# portion of the simulation box, and delta is the slip length determined by projecting the velocity profile of \n")
297 outFile.write("# liquid into the solid. \n#\n")
298 outFile.write("# zLiq = " + str(zLiq) + "\n")
299 outFile.write("# zSol = " + str(zSol) + "\n")
300 outFile.write("# delta = " + str(delta) + "\n")
301 if (jzType == 0):
302 JzPx = rnemdJzP[0]
303 outFile.write("# Jz(px) = " + str(JzPx) + "\n#\n")
304 elif (jzType == 1):
305 JzPy = rnemdJzP[1]
306 outFile.write("# Jz(py) = " + str(JzPy) + "\n#\n")
307 outFile.write("# shear rate \t\t\tlambda \n")
308 outFile.write("# " + str(shearRate) + "\t\t" + str(lambda1) + "\n#\n")
309 if (lambda1 < 0.0):
310 outFile.write("# WARNING NEGATIVE LAMBDA VALUE FOUND \n#\n#\n")
311
312
313def writeOutputKappa(outputFileName, rnemdJzP, jzType, zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2, popt, shearRate):
314 outFile = open(outputFileName, "a")
315 outFile.write("# The interfacial friction coefficient, kappa, was calculated by \n")
316 outFile.write("# \tkappa = Jz(p) / deltaV_{interface} \n")
317 outFile.write("# where Jz(p) is the imposed momentum flux, and deltaV_{interface} is the difference between the solid and \n")
318 outFile.write("# liquid velocities, obtained from the optimized fit measured across the interface \n#\n")
319 outFile.write("# \t left interface \t\t right interface \n")
320 outFile.write("# zLiq = " + str(zLiq1) + "\t\t\t" + str(zLiq2) + "\n")
321 outFile.write("# zSol = " + str(zSol1) + "\t\t\t" + str(zSol2) + "\n")
322
323 outFile.write("# vLiq = " + str(funcParabOnePoint(zLiq1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t\t" + str(funcParabOnePoint(zLiq2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]))+ "\n")
324
325 outFile.write("# vSol = " + str(funcParabOnePoint(zSol1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t\t" + str(funcParabOnePoint(zSol2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]))+ "\n")
326
327 outFile.write("# deltaV = " + str(deltaV1) + "\t\t" + str(deltaV2) + "\n")
328 ml = popt[3]
329 if (jzType == 0):
330 JzPx = rnemdJzP[0]
331 viscosity = 16.60539 * JzPx / ml
332 outFile.write("# Jz(px) = " + str(JzPx) + "\n#\n")
333 elif (jzType == 1):
334 JzPy = rnemdJzP[1]
335 viscosity = 16.60539 * JzPy / ml
336 outFile.write("# Jz(py) = " + str(JzPy) + "\n#\n")
337 outFile.write("# \t shear rate = " + str(shearRate) + "\n")
338 outFile.write("# \t kappa1 = " + str(kappa1) + "\n")
339 outFile.write("# \t kappa2 = " + str(kappa2) + "\n")
340 outFile.write("# \t viscosity = " + str(viscosity) + " mPa . s\n")
341 # outFile.write("# shear rate \t\t kappa1 \t\t kappa2 \n")
342 # outFile.write(str(shearRate) + "\t" + str(float(kappa1)) + "\t\t" + str(float(kappa2)) + "\n")
343 if (kappa1 < 0.0 or kappa2 < 0.0):
344 outFile.write("# WARNING NEGATIVE KAPPA VALUE FOUND \n")
345
346
347def main(argv):
348 parser = argparse.ArgumentParser(
349 description='OpenMD solid/liquid kinetic friction coefficient calculator for orthorhombic systems.',
350 #formatter_class=RawDescriptionHelpFormatter,
351 epilog="Example: solLiqFricCalc -i shearSim.rnemd -o shearSim.vfit -z1 30.4 -z2 75.6 -w 4.5 -Vs 1.2e-5 -Vl 1.0e-6 -m 2.2e-6 -d 2")
352 parser.add_argument("-i", "--rnemd-file=", action="store", dest="rnemdFileName", help="use specified input (.rnemd) file")
353 parser.add_argument("-o", "--vfit-file=", action="store", dest="vfitFileName", help="use specified output (.vfit) file")
354 parser.add_argument("-z1", "--lowerGibbsZ=", action="store", type=float, dest="z1", help="the location of the lower Gibbs dividing surface")
355 parser.add_argument("-z2", "--upperGibbsZ=", action="store", type=float, dest="z2", help="the location of the upper Gibbs dividing surface")
356 parser.add_argument("-l", "--lowerZVal=", action="store", nargs='?', type=float, dest="l", help="the initial estimate of the lower interface location (default=z1)")
357 parser.add_argument("-u", "--upperZVal=", action="store", nargs='?', type=float, dest="u", help="the initial estimate of the upper interface location (default=z2)")
358 parser.add_argument("-w", "--intWidth=", action="store", type=float, dest="w", help="the width of the interface")
359 parser.add_argument("-Vs", "--solidVel=", action="store", type=float, dest="Vs", help="the initial estimate of the velocity of the solid")
360 parser.add_argument("-Vl", "--liquidVel=", action="store", type=float, dest="Vl", help="the initial estimate of the velocity of the liquid")
361 parser.add_argument("-m", "--liquidSlope=", action="store", type=float, dest="m", help="the initial estimate of the slope in the liquid")
362 parser.add_argument("-d", "--toDelete=", action="store", default=0, type=int, dest="toDel", help="the number of data points to be deleted from the beginning and end of the velocity profile. (default=0)")
363 parser.add_argument("-s", "--sigma=", action="store", default=3.16549, type=float, dest="sigma", help="the molecular diameter of the liquid (default=3.16549)")
364 parser.add_argument("-f", "--convFactor=", action="store", default=2.19723, type=float, dest="convFactor", help="the conversion factor between widths obtained by fits and broader widths (default=2.19723 gives 90%%-10%% width)")
365 parser.add_argument("-t", "--useWidth=", action="store", default=False, type=bool, dest="useWidth", help="use provided width for deltaV calculation? (default=false)")
366 #parser.add_argument("-z","--useZLocations=", action="store", default=False, type=bool, dest="useZLocations", help="use provided z1 and z2 locations for calculations? (default=false)")
367
368 parser.add_argument("--useZLocations", action="store_true", dest="useZLocations", help="use the supplied z-locations z1 and z2 to calculate kappa")
369 parser.add_argument("--no-useZLocations", action="store_false", dest="useZLocations", help="use the fit values for z1 and z2 in the kappa calculation")
370 parser.set_defaults(useZLocations=True)
371
372 parser.add_argument("--lowerBound", action="store", type=float, dest="lowerBound", help="lower bound of the integral for Cf calculation")
373 parser.add_argument("--upperBound", action="store", type=float, dest="upperBound", help="upper bound of the integral for Cf calculation")
374 #parser.add_argument("--ms", action="store", type=float, dest="ms", help="slope of the velocity profile in the solid, should be close to zero.")
375 parser.add_argument("--boxlZ", action="store", dest="boxlZ", type=float, help="the box Z-dimension")
376 if len(sys.argv) == 1:
377 parser.print_help()
378 sys.exit(2)
379 args = parser.parse_args()
380
381 if (not args.rnemdFileName):
382 parser.print_help()
383 parser.error("No input-file was specified")
384
385 if (not args.vfitFileName):
386 parser.print_help()
387 parser.error("No output-file was specified")
388
389 if (not args.z1):
390 parser.print_help()
391 parser.error("No lower Gibbs dividing surface specified")
392
393 if (not args.z2):
394 parser.print_help()
395 parser.error("No upper Gibbs dividing surface specified")
396
397 if (not args.l):
398 args.l = args.z1
399
400 if (not args.u):
401 args.u = args.z2
402
403 if (not args.w):
404 parser.print_help()
405 parser.error("No width of interface specified")
406
407 if (not args.Vs):
408 parser.print_help()
409 parser.error("No initial estimate of the solid velocity specified")
410
411 if (not args.Vl):
412 parser.print_help()
413 parser.error("No initial estimate of the liquid velocity specified")
414
415 if (not args.m):
416 parser.print_help()
417 parser.error("No initial estimate of the slope in the liquid velocity profile specified")
418
419 #Call functions here, pass appropriate variables.
420
421 (rnemdJzP, jzType, rnemdZPos, rnemdV) = rnemdExtractor(args.rnemdFileName)
422 (popt, pcov, poptParab, pcovParab) = velFitter(rnemdJzP, jzType, rnemdZPos, rnemdV, args.Vs, args.l, args.u, args.m, args.w, args.Vl, args.toDel, args.vfitFileName, args.useZLocations)
423
424 (shearRate) = calcShearRate(jzType, rnemdZPos, rnemdV)
425 (zLiq, zSol, delta, lambda1) = lambdaCalc(rnemdJzP, jzType, poptParab, args.sigma, args.convFactor)
426 (zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2) = kappaCalc(rnemdJzP, jzType, args.z1, args.z2, args.w, poptParab, args.sigma, args.convFactor, args.useWidth, args.useZLocations)
427
428 writeOutputFile(args.vfitFileName, args.rnemdFileName, poptParab)
429 writeOutputLambda(args.vfitFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate)
430 writeOutputKappa(args.vfitFileName, rnemdJzP, jzType, zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2, poptParab, shearRate)
431
432if __name__ == "__main__":
433 main(sys.argv[1:])