OpenMD 3.1
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, least_squares
15from argparse import RawDescriptionHelpFormatter
16from scipy import stats
17
18def usage():
19 print(__doc__)
20
21# Trapezoidal sum for integrals
22def trapSum(A, B, fA, fB):
23 return 0.5*(B-A)*(fA+fB)
24
25# Fit function we are optimizing data to
26def func(x, Vs, Z1, Z2, ml, w, Vl):
27 ''' from Mathematica CForm
28 0.5*(Vl + ml*(x - Z1))*(1 - Tanh((x - Z1)/w)) +
29 0.5*Vs*(Tanh((x - Z1)/w) - Tanh((x - Z2)/w)) +
30 0.5*(Vl - ml*(x - Z2))*(1 + Tanh((x - Z2)/w))
31 '''
32
33 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)) )
34
35# Fit function we are optimizing data to
36def funcParab(x, Vs, Z1, Z2, ml, w, Vl):
37
38 y = np.empty_like (x)
39
40 k = 2.0 * (Vs - Vl - ml * (Z1-w)) / (w*w)
41
42 for i in range(len(x)):
43 xi = x[i]
44 if (xi < (Z1-w)):
45 yi = Vl + ml * xi
46 elif ((xi > Z1-w) and (xi < Z1) ):
47 yi = Vs - k * (xi - Z1)*(xi - Z1) / 2.0
48 elif ((xi > Z1) and (xi < Z2)):
49 yi = Vs
50 elif ((xi > Z2) and (xi < Z2 + w)):
51 yi = Vs - k * (xi - Z2)*(xi - Z2) / 2.0
52 else:
53 yi = Vs - k * w * w / 2.0 - ml * (xi - Z2 -w)
54 y[i] = yi
55
56 return y
57
58
59
60
61# first derivative of the fitting function
62def firstDeriv(x, Vs, Z1, Z2, ml, w, Vl):
63 ''' from Mathematica CForm
64 ((-0.5*Vl + 0.5*Vs - 0.5*ml*x + 0.5*ml*Z1)*Power(Sech((x - Z1)/w),2) +
65 (0.5*Vl - 0.5*Vs - 0.5*ml*x + 0.5*ml*Z2)*Power(Sech((x - Z2)/w),2) +
66 ml*w*(-0.5*Tanh((x - Z1)/w) - 0.5*Tanh((x - Z2)/w)))/w
67 '''
68
69 return ( ((-0.5*Vl + 0.5*Vs - 0.5*ml*x + 0.5*ml*Z1)*(1.0/np.cosh((x-Z1)/w))*(1.0/np.cosh((x-Z1)/w))+ (0.5*Vl - 0.5*Vs - 0.5*ml*x + 0.5*ml*Z2)*(1.0/np.cosh((x-Z2)/w))*(1.0/np.cosh((x-Z2)/w)) + ml*w*(-0.5*np.tanh((x - Z1)/w) - 0.5*np.tanh((x - Z2)/w)))/w )
70
71
72
73# Second derivative of the fit function
74def secondDeriv(x, Vs, Z1, Z2, ml, w, Vl):
75
76 ''' from Mathematica CForm
77 (Power(Sech((x - Z1)/w),2)*(-1.*m1*w +
78 (1.*Vl - 1.*Vs + 1.*m1*x - 1.*m1*Z1)*Tanh((x - Z1)/w)) +
79 Power(Sech((x - Z2)/w),2)*(-1.*m1*w +
80 (-1.*Vl + 1.*Vs + 1.*m1*x - 1.*m1*Z2)*Tanh((x - Z2)/w)))/Power(w,2)
81 '''
82
83 return ((1.0/np.cosh((x-Z1)/w))*(1.0/np.cosh((x-Z1)/w))*(-ml*w +
84 (Vl - Vs + ml*x - ml*Z1)*np.tanh((x - Z1)/w)) +
85 (1.0/np.cosh((x-Z2)/w))*(1.0/np.cosh((x-Z2)/w))*(-ml*w +
86 (-Vl + Vs + ml*x - ml*Z2)*np.tanh((x - Z2)/w)))/(w*w)
87
88
89
90
91
92
93
94
95
96
97#This function reads in the .rnemd file and extracts the actual momentum flux, as well as the zPositions and the average velocities therein.
98def rnemdExtractor(rnemdFileName):
99 if os.path.exists(rnemdFileName):
100 rnemdFile = open(rnemdFileName, 'r')
101 else:
102 #print 'Error: cannot open ' + rnemdFileName
103 sys.exit(2)
104
105 rnemdZPos = np.array([])
106 rnemdV = np.array([]).reshape(0, 3)
107
108 while True:
109 line = rnemdFile.readline()
110 if not line: break
111 if "Actual flux:" in line or "actual flux:" in line:
112 rnemdFile.readline()
113 line = rnemdFile.readline()
114 rnemdJzP = np.array([[float(line.split()[4][:-1])], [float(line.split()[5][:-1])], [float(line.split()[6])]])
115 elif "#" not in line:
116 if len(line.split()) > 6:
117 print("Error in the .rnemd file, found more than 6 expected columns.")
118 else:
119 rnemdZPos = np.append(rnemdZPos, np.array([float(line.split()[0])]))
120 rnemdV = np.vstack( [rnemdV, [ float(line.split()[2]), float(line.split()[3]), float(line.split()[4]) ] ] )
121
122 #Now let's determine what kind of shearing simulation was performed
123 # ie, JzPx, JzPy, or JzPxy
124 if (rnemdJzP[2] != 0.0):
125 print("Error: JzPz Actual Exchange total was found to be >0.0 in (.rnemd) file.")
126 sys.exit(2)
127 if (rnemdJzP[0] == 0.0 and rnemdJzP[1] == 0.0):
128 print("Error: both JzPx and JzPy Actual Exchange totals were found to be 0.0 in (.rnemd) file. ")
129 sys.exit(2)
130 elif (rnemdJzP[0] != 0.0 and rnemdJzP[1] != 0.0):
131 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.")
132 sys.exit(2)
133 #After sanity checking, we now set a variable jzType to store if the imposed momentum flux is in the x- (0) or y-direction (1). #jzType will be used when we iterate through rnemdV
134 # jzType = 0 -> flux was in x-dimension
135 # jzType = 1 -> flux was in y-dimension
136 elif (rnemdJzP[0] == 0.0 or rnemdJzP[1] == 0.0):
137 if np.abs(rnemdJzP[0]) > np.abs(rnemdJzP [1]):
138 jzType = 0
139 else:
140 jzType = 1
141
142 #print "Read in " + str(rnemdFileName)
143 return (rnemdJzP, jzType, rnemdZPos, rnemdV)
144
145
146
147
148
149
150
151
152def velFitter(rnemdJzP, jzType, rnemdZPos, rnemdV, Vs, Z1, Z2, m, w, Vl, numDel, outputFileName, useZLocations):
153 rnemdVx = np.array([]).reshape(0, 1)
154 rnemdVy = np.array([]).reshape(0, 1)
155
156 for i in range(0, len(rnemdV)):
157 rnemdVx = np.append(rnemdVx, rnemdV[i][0])
158 rnemdVy = np.append(rnemdVy, rnemdV[i][1])
159
160 #Remove the first two and last two values of rnemdVx and rnemdVy to smooth fits
161 for i in range(0, numDel):
162 rnemdVx = np.delete(rnemdVx, 0)
163 rnemdVy = np.delete(rnemdVy, 0)
164 rnemdZPos = np.delete(rnemdZPos, 0)
165 for i in range(0, numDel+1):
166 rnemdVx = np.delete(rnemdVx, len(rnemdVx)-1)
167 rnemdVy = np.delete(rnemdVy, len(rnemdVy)-1)
168 rnemdZPos = np.delete(rnemdZPos, len(rnemdZPos)-1)
169
170 #Tanh fitting function...
171 #curve_fit(function to fit to, x-data, y-data, [initial guess of param 1, initial guess of param 2,...]
172 if (jzType == 0):
173 popt, pcov = least_squares(func, rnemdZPos, rnemdVx, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
174 elif (jzType == 1):
175 popt, pcov = curve_fit(func, rnemdZPos, rnemdVy, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
176
177
178 #write out the raw data and the fit for visual inspection of quality of fit
179 #print "Writing out fit file"
180 fitFileName = outputFileName + "FIT"
181 fitFile = open(fitFileName, "w")
182 fitFile.write("# zPosition <V_{i}> fit value" + "\n")
183 for i in range(0, len(rnemdZPos)):
184 y = func(rnemdZPos[i], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
185 if (jzType == 0):
186 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVx[i]) + "\t" + str(y) + "\n")
187 elif (jzType == 1):
188 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVy[i]) + "\t" + str(y) + "\n")
189
190 # Parabola fitting function...
191 #curve_fit(function to fit to, x-data, y-data, [initial guess of param 1, initial guess of param 2,...]
192 # funcParab(x, Vs, Z1, Z2, ml, w, Vl)
193 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])
194 if (jzType == 0):
195 poptParab, pcovParab = curve_fit(funcParab, rnemdZPos, rnemdVx, [Vs, Z1, Z2, m, w, Vl], bounds=param_bounds, ftol=0.01)
196 elif (jzType == 1):
197 poptParab, pcovParab = curve_fit(funcParab, rnemdZPos, rnemdVy, [Vs, Z1, Z2, m, w, Vl], bounds=param_bounds, ftol=0.01)
198
199
200
201 #print "popt = ", popt
202 return (popt, pcov, poptParab, pcovParab)
203
204
205
206
207
208
209
210def calcShearRate(jzType, rnemdZPos, rnemdV):
211 # jzType = 0 -> flux was in x-dimension # jzType = 1 -> flux was in y-dimension
212
213 velLiq = (rnemdV[0][jzType] + rnemdV[len(rnemdV)-1][jzType])/2.0
214
215 if (len(rnemdV) % 2) == 0 :
216 velSol = ( rnemdV[(len(rnemdV) / 2)-1][jzType] + rnemdV[(len(rnemdV)/ 2)][jzType]) / 2.0
217 elif (len(rnemdV) % 2) != 0:
218 velSol = ( rnemdV[(len(rnemdV) + 1)/2][jzType] + rnemdV[((len(rnemdV)+1)/2)+1][jzType]) / 2.0
219
220 shearRate = velSol - velLiq
221
222 return (shearRate)
223
224
225
226
227
228
229def lambdaCalc(rnemdJzP, jzType, popt, sigma, convFactor):
230 #print "Performing lambda calculation"
231 #popt = Vs, Z1, Z2, ml, w, Vl
232
233 zLiq = (popt[0] - popt[5]) / popt[3] + popt[1]
234 if (popt[4]*convFactor < sigma):
235 zSol = popt[1] + 0.5*sigma
236 else:
237 zSol = popt[1] + 0.5*(popt[4]*convFactor)
238 delta = zLiq - zSol
239
240 if (jzType == 0):
241 lambda1 = rnemdJzP[0] / (delta * popt[3])
242 elif (jzType == 1):
243 lambda1 = rnemdJzP[1] / (delta * popt[3])
244
245
246 return (zLiq, zSol, delta, lambda1)
247
248
249
250
251
252
253def kappaCalc(rnemdJzP, jzType, z1, z2, w, popt, sigma, convFactor, useWidth, useZLocations):
254 #print "Performing kappa calculation"
255
256 #print "useZLocations = " , useZLocations
257
258 if (useZLocations):
259 print("Using supplied locations of the interface for kappa")
260 if (useWidth):
261 zLiq1 = z1 - 0.5*w
262 zSol1 = z1 + 0.5*w
263 zLiq2 = z2 + 0.5*w
264 zSol2 = z2 - 0.5*w
265 elif (not useWidth):
266 zLiq1 = z1 - 1.5*sigma
267 zSol1 = z1 + 1.5*sigma
268 zLiq2 = z2 + 1.5*sigma
269 zSol2 = z2 - 1.5*sigma
270 elif (not useZLocations):
271 print("Using the fit values for the locations of the interface for kappa")
272 if (useWidth):
273 zLiq1 = popt[1] - 0.5*w
274 zSol1 = popt[1] + 0.5*w
275 zLiq2 = popt[2] + 0.5*w
276 zSol2 = popt[2] - 0.5*w
277 elif (not useWidth):
278 zLiq1 = popt[1] - 1.5*sigma
279 zSol1 = popt[1] + 1.5*sigma
280 zLiq2 = popt[2] + 1.5*sigma
281 zSol2 = popt[2] - 1.5*sigma
282
283 vLiq1 = func(zLiq1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
284 vSol1 = func(zSol1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
285 deltaV1 = vSol1 - vLiq1
286
287 vLiq2 = func(zLiq2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
288 vSol2 = func(zSol2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
289 deltaV2 = vSol2 - vLiq2
290
291
292 if (jzType == 0):
293 k1 = rnemdJzP[0] / deltaV1
294 k2 = rnemdJzP[0] / deltaV2
295 elif (jzType == 1):
296 k1 = rnemdJzP[1] / deltaV1
297 k2 = rnemdJzP[1] / deltaV2
298
299 return (zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, k1, k2)
300
301
302
303
304
305
306
307def secondDerivVelApprox(rnemdJzP, jzType, z1, z2, w, popt, useZLocations):
308 #print "Performing second derivative kappa calculation"
309
310 if (useZLocations):
311 zLiq1 = z1 - 0.5*w
312 zSol1 = z1 + 0.5*w
313 vInt1 = func(z1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
314
315 zLiq2 = z2 + 0.5*w
316 zSol2 = z2 - 0.5*w
317 vInt2 = func(z2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
318
319 elif (not useZLocations):
320 zLiq1 = popt[1] - 0.5*w
321 zSol1 = popt[1] + 0.5*w
322 vInt1 = func(popt[1], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
323
324 zLiq2 = popt[2] + 0.5*w
325 zSol2 = popt[2] - 0.5*w
326 vInt2 = func(popt[2], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
327
328
329 # using the fit z1 and z2 values, the kappa's should be the same
330 # due to the symmetric nature of the fit.
331 vLiq1 = func(zLiq1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
332 vSol1 = func(zSol1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
333
334 vLiq2 = func(zLiq2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
335 vSol2 = func(zSol2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
336
337
338 #Testing parameters
339 jzType = 0
340 rnemdJzP[0] = 1.0
341 vSol1 = 1.0
342 vInt1 = 1.0
343 vLiq1 = 1.0
344
345 if (jzType == 0):
346 momentumConduct1 = (4.0*rnemdJzP[0]) * (vSol1 - 2*vInt1 + vLiq1) / (vSol1*vSol1 - 2*vSol1*vLiq1 + vLiq1*vLiq1)
347 momentumConduct2 = (4.0*rnemdJzP[0]) * (vSol2 - 2*vInt2 + vLiq2) / (vSol2*vSol2 - 2*vSol2*vLiq2 + vLiq2*vLiq2)
348
349 elif (jzType == 1):
350 momentumConduct1 = (4.0*rnemdJzP[1]) * (vSol1 - 2*vInt1 + vLiq1) / (vSol1*vSol1 - 2*vSol1*vLiq1 + vLiq1*vLiq1)
351 momentumConduct2 = (4.0*rnemdJzP[1]) * (vSol2 - 2*vInt2 + vLiq2) / (vSol2*vSol2 - 2*vSol2*vLiq2 + vLiq2*vLiq2)
352
353 #print "momentumConduct1 = " , momentumConduct1
354 #print "momentumConduct2 = " , momentumConduct2
355
356
357
358
359
360
361
362
363def writeOutputFile(outputFileName, rnemdFileName, popt):
364 outFile = open(outputFileName, "w")
365
366 outFile.write("##################################################### \n")
367 outFile.write("## This output file was generated by velocityFitter.py on " + str(datetime.datetime.now()) + " ## \n#\n")
368 outFile.write("# The velocity profile found in " + str(rnemdFileName) + " was fit by \n")
369 outFile.write("# \ty = (Vl + ml(x-Z1))*(0.5(1-tanh((x-Z1)/w))) + Vs*0.5(tanh((x-Z1)/2)-tanh((x-Z2)/w)) + (Vl-ml(x-Z2))*(0.5(1+tanh((x-Z2/w)))) \n")
370 outFile.write("# where (Vs, Vl, Z1, Z2, ml, w) are fit parameters, and here y = V(i) and x = z \n")
371 outFile.write("# Vs = velocity of the solid\n# Vl = velocity of the liquid at the interface \n")
372 outFile.write("# Z1 = z-position of the lower interface \n# Z2 = z-position of the upper interface \n")
373 outFile.write("# ml = slope of the velocity profile in the lower liquid region of the box \n")
374
375 outFile.write("# w = width of the interface \n#\n")
376 outFile.write("# Obtained optomized parameters \n")
377 outFile.write("# \t Vs = " + str(popt[0]) + "\n")
378 outFile.write("# \t Vl = " + str(popt[5]) + "\n")
379 outFile.write("# \t Z1 = " + str(popt[1]) + "\n")
380 outFile.write("# \t Z2 = " + str(popt[2]) + "\n")
381 outFile.write("# \t ml = " + str(popt[3]) + "\n")
382 outFile.write("# \t w = " + str(popt[4]) + "\n#\n")
383
384 if (float(popt[0]) < float(popt[5])):
385 print("Bad fit detected, Vs < Vl")
386 outFile.write("# WARNING Vs < Vl \n#\n")
387
388
389
390
391
392def writeOutputLambda(outputFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate):
393 outFile = open(outputFileName, "a")
394 outFile.write("# The interfacial friction coefficient, lambda, was calculated by \n")
395 outFile.write("# \tlambda = Jz(p) / (dV/dz)_{liquid} / delta \n")
396 outFile.write("# where Jz(p) is the imposed momentum flux, (dV/dz)_{liquid} is the slope of the velocity profile in the liquid \n")
397 outFile.write("# portion of the simulation box, and delta is the slip length determined by projecting the velocity profile of \n")
398 outFile.write("# liquid into the solid. \n#\n")
399 outFile.write("# zLiq = " + str(zLiq) + "\n")
400 outFile.write("# zSol = " + str(zSol) + "\n")
401 outFile.write("# delta = " + str(delta) + "\n")
402 if (jzType == 0):
403 outFile.write("# Jz(px) = " + str((float(rnemdJzP[0]))) + "\n#\n")
404 elif (jzType == 1):
405 outFile.write("# Jz(py) = " + str((float(rnemdJzP[1]))) + "\n#\n")
406 outFile.write("# shear rate \t\t\tlambda \n")
407 outFile.write("# " + str(shearRate) + "\t\t" + str(float(lambda1)) + "\n#\n")
408 if (lambda1 < 0.0):
409 outFile.write("# WARNING NEGATIVE LAMBDA VALUE FOUND \n#\n#\n")
410
411
412
413
414
415
416def writeOutputKappa(outputFileName, rnemdJzP, jzType, zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2, popt, shearRate):
417 outFile = open(outputFileName, "a")
418 outFile.write("# The interfacial friction coefficient, kappa, was calculated by \n")
419 outFile.write("# \tkappa = Jz(p) / deltaV_{interface} \n")
420 outFile.write("# where Jz(p) is the imposed momentum flux, and deltaV_{interface} is the difference between the solid and \n")
421 outFile.write("# liquid velocities, obtained from the optimized fit measured across the interface \n#\n")
422 outFile.write("# \t left interface \t\t right interface \n")
423 outFile.write("# zLiq = " + str(zLiq1) + "\t\t\t" + str(zLiq2) + "\n")
424 outFile.write("# zSol = " + str(zSol1) + "\t\t\t" + str(zSol2) + "\n")
425
426 outFile.write("# vLiq = " + str(func(zLiq1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t\t" + str(func(zLiq2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]))+ "\n")
427
428 outFile.write("# vSol = " + str(func(zSol1, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t\t" + str(func(zSol2, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]))+ "\n")
429
430 outFile.write("# deltaV = " + str(deltaV1) + "\t\t" + str(deltaV2) + "\n")
431 if (jzType == 0):
432 outFile.write("# Jz(px) = " + str((float(rnemdJzP[0]))) + "\n#\n")
433 elif (jzType == 1):
434 outFile.write("# Jz(py) = " + str((float(rnemdJzP[1]))) + "\n#\n")
435 outFile.write("# \t shear rate = " + str(shearRate) + "\n")
436 outFile.write("# \t kappa1 = " + str(kappa1) + "\n")
437 outFile.write("# \t kappa2 = " + str(kappa2) + "\n")
438 #outFile.write("# shear rate \t\t kappa1 \t\t kappa2 \n")
439 #outFile.write(str(shearRate) + "\t" + str(float(kappa1)) + "\t\t" + str(float(kappa2)) + "\n")
440 if (kappa1 < 0.0 or kappa2 < 0.0):
441 outFile.write("# WARNING NEGATIVE KAPPA VALUE FOUND \n")
442
443
444
445def writeOutputSecondDeriv(outputFileName, popt):
446 outFile = open(outputFileName, "a")
447
448 outFile.write("#\n#\n#")
449 outFile.write("# second derivative of fit at Z1 \n")
450 outFile.write("# Z1 = \t " + str(popt[1]) + "\n")
451 outFile.write("# d^2fit/dz^2 |z1 \n")
452 outFile.write(str(secondDeriv(popt[1], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
453
454
455
456
457
458def writeOutputForceMethod(outputFileName, rnemdJzP, jzType, popt):
459 outFile = open(outputFileName, "a")
460
461
462 outFile.write("#\n#\n")
463 outFile.write("# The Force matching method is reported here. We assume the system has come to steady state at which point the Fappl = sum(Ffric) \n")
464 outFile.write("# \t Fsl = jz(p) + \eta * (dV/dz) \n")
465 outFile.write("# \t Cf = 0.5*[1 + \eta * (dV/dz) / jz(p)] \n")
466
467 viscosity = 0.88 # Jz(px) / (dVx/dz) from bulk SPC/E @ 225K
468
469 if (jzType == 0):
470 Fsl = rnemdJzP[0] + viscosity*popt[3]
471 Cf = 0.5 + (0.5*viscosity*popt[3]/rnemdJzP[0])
472 outFile.write("# jz(px) = " + str(float(rnemdJzP[0])) + "\n")
473 outFile.write("# \eta = " + str(viscosity) + "\n")
474 outFile.write("# (dVx/dz) = " + str(popt[3]) + "\n")
475
476 elif (jzType == 1):
477 Fsl = rnemdJzP[1] + viscosity*popt[3]
478 Cf = 0.5 + (0.5*viscosity*popt[3]/rnemdJzP[1])
479 outFile.write("# jz(py) = " + str(float(rnemdJzP[1])) + "\n")
480 outFile.write("# \eta = " + str(viscosity) + "\n")
481 outFile.write("# (dVy/dz) = " + str(popt[3]) + "\n#\n")
482
483
484 outFile.write("# Fsl = " + str(float(Fsl)) + "\n")
485 outFile.write("# Cf_{ForceMethod} = " + str(float(Cf)) + "\n")
486
487
488
489
490
491def writeOutputIntegralMethod(outputFileName, rnemdJzP, jzType, popt, lowerBound, upperBound, boxlZ, poptParab):
492 derivFileName = outputFileName + "DERIV"
493 derivFile = open(derivFileName, "w")
494 derivFile.write("# xPos \t U \t dU/dz \t d^2U/dz^2 \t (d^2U/dz^2)\‍(dU/dz) \n")
495
496 outFile = open(outputFileName, "a")
497
498 outFile.write("#\n#\n")
499 outFile.write("# The integral method makes no assumptions about the system having reached \n")
500 outFile.write("# a steady-state. \n")
501 outFile.write("# \t Cf = F_{fric} \ F_{appl} \n")
502 outFile.write("# \t = 2* Integral[ (d^2u/dz^2) / (du/dz) dz]|^{solid}_{liquid} \n")
503 outFile.write("# Integral bounds: \n")
504 outFile.write("# \t lowerBound = " + str(lowerBound) + "\n")
505 outFile.write("# \t upperBound = " + str(upperBound) + "\n")
506
507
508
509 # print out the derivatives for visual inspection
510 deltaX = boxlZ/1000
511 xPos = 0.0
512 while (xPos < boxlZ):
513 derivFile.write(str(xPos) + "\t" + str(func(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t" + str(firstDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t" + str(secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\t" + str(secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])/firstDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
514 xPos = xPos + deltaX
515
516
517 '''
518 # Write out smooth porabola fit
519 x = np.array([])
520 xPos = 0.0
521 while (xPos < boxlZ):
522 x = np.append(x, np.array([xPos]))
523 xPos = xPos + deltaX
524
525 parabFileName = outputFileName + "Parab"
526 parabFile = open(parabFileName, "w")
527 parabY = funcParab(x, poptParab[0], poptParab[1], poptParab[2], poptParab[3], poptParab[4], poptParab[5])
528
529 for i in range(0, len(x)):
530 parabFile.write(str(x[i]) + "\t" + str(parabY[i]) + "\n")
531
532 #def funcParab(x, Vs, Z1, Z2, ml, w, Vl):
533 #k = 2.0 * (Vs - Vl - ml * (Z1-w)) / (w*w)
534 outFile.write("# Parabola Fitting Method \n")
535 outFile.write("# Cf = 4kw\mu \n#\n")
536 outFile.write("# Parameters: \n")
537 k = 2.0 * (poptParab[0] - poptParab[5] - poptParab[3] * (poptParab[1] - poptParab[4])) / ( poptParab[4] * poptParab[4])
538 outFile.write("# k = " + str(k) + "\n")
539 outFile.write("# w = " + str(poptParab[4]) + "\n")
540 outFile.write("# \mu = " + str(1) + "\n")
541 CfParab = 4.0 * k * 1.0 * poptParab[4]
542 outFile.write("# Cf_{Parab} = " + str(CfParab) + "\n")
543 '''
544
545
546 # calculate integral Cf
547 deltaX = (upperBound - lowerBound)/1000
548 integral = 0.0
549 xPos = lowerBound
550 while (xPos < upperBound):
551 lowerTerm = secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]) / firstDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
552 upperTerm = secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]) / firstDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
553 dArea = trapSum(xPos, xPos + deltaX, lowerTerm, upperTerm)
554 integral = integral + dArea
555 xPos = xPos + deltaX
556
557 integral = -2.0 * integral #due to the 2 interfaces present in systems
558
559 outFile.write("# Cf_{IntegralMethod} = " + str(integral) + "\n")
560
561
562
563
564
565def writeOutputSimpleCfCalc(outputFileName, rnemdJzP, jzType, popt, givenWidth, boxlZ, poptParab):
566
567 # Find the minimum of (d^2U/dz^2)
568 xPos = 0.0
569 deltaX = boxlZ/1000
570 minUpp = 1e6
571 minUppXPos = 0.0
572 while (xPos < (boxlZ/2.0)):
573 if (secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
574 < minUpp):
575 minUpp = secondDeriv(xPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
576 minUppXPos = xPos
577 xPos = xPos + deltaX
578
579
580 outFile = open(outputFileName, "a")
581 outFile.write("#\n#\n#\n")
582 outFile.write("# Simple Cf Calculation with fit width \n")
583 outFile.write("# \t Cf = -2.0*w (d^2u\dz^2)|z=Z1 / (du\dz)|z=Z1 \n")
584 outFile.write("#\n# Parameters in calculation \n")
585 outFile.write("# w = " + str(popt[4]) + "\n")
586 outFile.write("# (d^2u\dz^2)|z=Z1 = " + str(secondDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
587 outFile.write("# (du\dz)|z=Z1 = " + str(firstDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
588
589 Cf = -2.0*popt[4]*secondDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])/ firstDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
590
591 outFile.write("#\n# Cf_{SimpleMethod_FitWidth} = " + str(Cf) + "\n")
592
593 '''
594 # Find the minimum of (d^2U/dz^2)
595 xPos = 0.0
596 deltaX = boxlZ/1000
597 minUpp = 1e6
598 minUppXPos = 0.0
599 while (xPos < (boxlZ/2.0)):
600 if (secondDeriv(xPos,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5])
601 < minUpp):
602 minUpp = secondDeriv(xPos,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5])
603 minUppXPos = xPos
604 xPos = xPos + deltaX
605 '''
606
607 # let's also calculate Cf_{SimpleMethod} with the width given by the user
608 outFile.write("#\n#\n#\n")
609 outFile.write("# Simple Cf Calculation with given width \n")
610 outFile.write("# \t Cf = -2.0*w (d^2u\dz^2)|z_{min(U'')} / (du\dz)|z_{min(U'')} \n")
611 outFile.write("#\n# Parameters in calculation \n")
612 outFile.write("# w = " + str(givenWidth) + "\n")
613 outFile.write("# minUppxPos = " + str(minUppXPos) + "\n")
614 outFile.write("# (d^2u\dz^2)|z_{min(U'')} = " + str(secondDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
615 outFile.write("# (du\dz)|z_{min(U'')} = " + str(firstDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) +"\n")
616 outFile.write("# (d^2u/dz^2)/(du/dz) = " + str(secondDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])/firstDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])) + "\n")
617
618 Cf = -2.0*givenWidth*secondDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])/ firstDeriv(minUppXPos, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
619
620 outFile.write("#\n# Cf_{SimpleMethod_GivenWidth} = " + str(Cf) + "\n")
621
622
623def ParabolaCf(outputFileName, rnemdJzP, jzType, popt, lowerBound, upperBound, boxlZ, poptParab):
624
625
626 # Write out smooth porabola fit
627 x = np.array([])
628 deltaX = boxlZ/1000
629 xPos = 0.0
630 while (xPos < boxlZ):
631 x = np.append(x, np.array([xPos]))
632 xPos = xPos + deltaX
633
634 parabFileName = outputFileName + "Parab"
635 parabFile = open(parabFileName, "w")
636 parabY = funcParab(x, poptParab[0], poptParab[1], poptParab[2], poptParab[3], poptParab[4], poptParab[5])
637
638 for i in range(0, len(x)):
639 parabFile.write(str(x[i]) + "\t" + str(parabY[i]) + "\n")
640
641
642 outFile = open(outputFileName, "a")
643 outFile.write("#\n#\n#\n")
644
645 #def funcParab(x, Vs, Z1, Z2, ml, w, Vl):
646 #k = 2.0 * (Vs - Vl - ml * (Z1-w)) / (w*w)
647 outFile.write("# Parabola Fitting Method \n")
648 outFile.write("# Cf = 2 * k * w / ml \n#\n")
649 outFile.write("# Parameters: \n")
650 k = 2.0 * (poptParab[0] - poptParab[5] - poptParab[3] * (poptParab[1] - poptParab[4])) / ( poptParab[4] * poptParab[4])
651 outFile.write("# k = " + str(k) + "\n")
652 outFile.write("# w = " + str(poptParab[4]) + "\n")
653 outFile.write("# ml = " + str(poptParab[3]) + "\n")
654 CfParab = 2.0 * k * poptParab[4] / poptParab[3]
655 outFile.write("# Cf_{Parab} = " + str(CfParab) + "\n")
656
657
658
659
660def main(argv):
661 parser = argparse.ArgumentParser(
662 description='OpenMD solid/liquid kinetic friction coefficient calculator for orthorhombic systems.',
663 #formatter_class=RawDescriptionHelpFormatter,
664 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")
665 parser.add_argument("-i", "--rnemd-file=", action="store", dest="rnemdFileName", help="use specified input (.rnemd) file")
666 parser.add_argument("-o", "--vfit-file=", action="store", dest="vfitFileName", help="use specified output (.vfit) file")
667 parser.add_argument("-z1", "--lowerGibbsZ=", action="store", type=float, dest="z1", help="the location of the lower Gibbs dividing surface")
668 parser.add_argument("-z2", "--upperGibbsZ=", action="store", type=float, dest="z2", help="the location of the upper Gibbs dividing surface")
669 parser.add_argument("-l", "--lowerZVal=", action="store", nargs='?', type=float, dest="l", help="the initial estimate of the lower interface location (default=z1)")
670 parser.add_argument("-u", "--upperZVal=", action="store", nargs='?', type=float, dest="u", help="the initial estimate of the upper interface location (default=z2)")
671 parser.add_argument("-w", "--intWidth=", action="store", type=float, dest="w", help="the width of the interface")
672 parser.add_argument("-Vs", "--solidVel=", action="store", type=float, dest="Vs", help="the initial estimate of the velocity of the solid")
673 parser.add_argument("-Vl", "--liquidVel=", action="store", type=float, dest="Vl", help="the initial estimate of the velocity of the liquid")
674 parser.add_argument("-m", "--liquidSlope=", action="store", type=float, dest="m", help="the initial estimate of the slope in the liquid")
675 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)")
676 parser.add_argument("-s", "--sigma=", action="store", default=3.16549, type=float, dest="sigma", help="the molecular diameter of the liquid (default=3.16549)")
677 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)")
678 parser.add_argument("-t", "--useWidth=", action="store", default=False, type=bool, dest="useWidth", help="use provided width for deltaV calculation? (default=false)")
679 #parser.add_argument("-z","--useZLocations=", action="store", default=False, type=bool, dest="useZLocations", help="use provided z1 and z2 locations for calculations? (default=false)")
680
681 parser.add_argument("--useZLocations", action="store_true", dest="useZLocations", help="use the supplied z-locations z1 and z2 to calculate kappa")
682 parser.add_argument("--no-useZLocations", action="store_false", dest="useZLocations", help="use the fit values for z1 and z2 in the kappa calculation")
683 parser.set_defaults(useZLocations=True)
684
685 parser.add_argument("--lowerBound", action="store", type=float, dest="lowerBound", help="lower bound of the integral for Cf calculation")
686 parser.add_argument("--upperBound", action="store", type=float, dest="upperBound", help="upper bound of the integral for Cf calculation")
687 #parser.add_argument("--ms", action="store", type=float, dest="ms", help="slope of the velocity profile in the solid, should be close to zero.")
688 parser.add_argument("--boxlZ", action="store", dest="boxlZ", type=float, help="the box Z-dimension")
689 if len(sys.argv) == 1:
690 parser.print_help()
691 sys.exit(2)
692 args = parser.parse_args()
693
694 if (not args.rnemdFileName):
695 parser.print_help()
696 parser.error("No input-file was specified")
697
698 if (not args.vfitFileName):
699 parser.print_help()
700 parser.error("No output-file was specified")
701
702 if (not args.z1):
703 parser.print_help()
704 parser.error("No lower Gibbs dividing surface specified")
705
706 if (not args.z2):
707 parser.print_help()
708 parser.error("No upper Gibbs dividing surface specified")
709
710 if (not args.l):
711 args.l = args.z1
712
713 if (not args.u):
714 args.u = args.z2
715
716 if (not args.w):
717 parser.print_help()
718 parser.error("No width of interface specified")
719
720 if (not args.Vs):
721 parser.print_help()
722 parser.error("No initial estimate of the solid velocity specified")
723
724 if (not args.Vl):
725 parser.print_help()
726 parser.error("No initial estimate of the liquid velocity specified")
727
728 if (not args.m):
729 parser.print_help()
730 parser.error("No initial estimate of the slope in the liquid velocity profile specified")
731
732
733
734
735
736
737 #Call functions here, pass appropriate variables.
738
739
740 (rnemdJzP, jzType, rnemdZPos, rnemdV) = rnemdExtractor(args.rnemdFileName)
741 (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)
742
743
744
745 (shearRate) = calcShearRate(jzType, rnemdZPos, rnemdV)
746 (zLiq, zSol, delta, lambda1) = lambdaCalc(rnemdJzP, jzType, popt, args.sigma, args.convFactor)
747 (zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2) = kappaCalc(rnemdJzP, jzType, args.z1, args.z2, args.w, popt, args.sigma, args.convFactor, args.useWidth, args.useZLocations)
748
749
750
751 #secondDerivVelApprox(rnemdJzP, jzType, args.z1, args.z2, args.w, popt, args.useZLocations)
752
753
754
755 writeOutputFile(args.vfitFileName, args.rnemdFileName, popt)
756 writeOutputLambda(args.vfitFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate)
757 writeOutputKappa(args.vfitFileName, rnemdJzP, jzType, zLiq1, zSol1, deltaV1, zLiq2, zSol2, deltaV2, kappa1, kappa2, popt, shearRate)
758
759 #def secondDeriv(x, Vs, Z1, Z2, ml, w, Vl):
760 # need to pass popt[1] = Z1 = x to obtain secondDeriv at the Z1 location
761 #writeOutputSecondDeriv(args.vfitFileName, popt)
762
763 #def calcForceMethod()
764 writeOutputForceMethod(args.vfitFileName, rnemdJzP, jzType, popt)
765
766
767 writeOutputIntegralMethod(args.vfitFileName, rnemdJzP, jzType, popt, args.lowerBound, args.upperBound, args.boxlZ, poptParab)
768 writeOutputSimpleCfCalc(args.vfitFileName, rnemdJzP, jzType, popt, args.w, args.boxlZ, poptParab)
769
770 ParabolaCf(args.vfitFileName, rnemdJzP, jzType, popt, args.lowerBound, args.upperBound, args.boxlZ, poptParab)
771
772
773
774if __name__ == "__main__":
775 main(sys.argv[1:])