OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
slipLength
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
18
19
20def usage():
21 print(__doc__)
22
23
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
36
37
38
39
40#This function reads in the .rnemd file and extracts the actual momentum flux, as well as the zPositions and the average velocities therein.
41def rnemdExtractor(rnemdFileName):
42 if os.path.exists(rnemdFileName):
43 rnemdFile = open(rnemdFileName, 'r')
44 else:
45 print('Error: cannot open ' + rnemdFileName)
46 sys.exit(2)
47
48 rnemdZPos = np.array([])
49 rnemdV = np.array([]).reshape(0, 3)
50
51 while True:
52 line = rnemdFile.readline()
53 if not line: break
54 if "Actual flux:" in line or "actual flux:" in line:
55 rnemdFile.readline()
56 line = rnemdFile.readline()
57 rnemdJzP = np.array([[float(line.split()[4][:-1])], [float(line.split()[5][:-1])], [float(line.split()[6])]])
58 elif "#" not in line:
59 if len(line.split()) > 6:
60 print("Error in the .rnemd file, found more than 6 expected columns.")
61 else:
62 rnemdZPos = np.append(rnemdZPos, np.array([float(line.split()[0])]))
63 rnemdV = np.vstack( [rnemdV, [ float(line.split()[2]), float(line.split()[3]), float(line.split()[4]) ] ] )
64
65 #Now let's determine what kind of shearing simulation was performed
66 # ie, JzPx, JzPy, or JzPxy
67 if (rnemdJzP[2] != 0.0):
68 print("Error: JzPz Actual Exchange total was found to be >0.0 in (.rnemd) file.")
69 sys.exit(2)
70 if (rnemdJzP[0] == 0.0 and rnemdJzP[1] == 0.0):
71 print("Error: both JzPx and JzPy Actual Exchange totals were found to be 0.0 in (.rnemd) file. ")
72 sys.exit(2)
73 elif (rnemdJzP[0] != 0.0 and rnemdJzP[1] != 0.0):
74 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.")
75 sys.exit(2)
76
77 #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).
78 #jzType will be used when we iterate through rnemdV
79 #jzType = 0 -> flux was in x-dimension
80 #jzType = 1 -> flux was in y-dimension
81 elif (rnemdJzP[0] == 0.0 or rnemdJzP[1] == 0.0):
82 if np.abs(rnemdJzP[0]) > np.abs(rnemdJzP [1]):
83 jzType = 0
84 else:
85 jzType = 1
86
87 print("Read in " + str(rnemdFileName))
88 return (rnemdJzP, jzType, rnemdZPos, rnemdV)
89
90
91
92
93
94
95
96
97def velFitter(rnemdJzP, jzType, rnemdZPos, rnemdV, Vs, Z1, Z2, m, w, Vl, numDel, outputFileName):
98 rnemdVx = np.array([]).reshape(0, 1)
99 rnemdVy = np.array([]).reshape(0, 1)
100
101 for i in range(0, len(rnemdV)):
102 rnemdVx = np.append(rnemdVx, rnemdV[i][0])
103 rnemdVy = np.append(rnemdVy, rnemdV[i][1])
104
105 #Remove data points from the start and last end of the vector of values of rnemdVx and rnemdVy to smooth fits
106 for i in range(0, numDel):
107 rnemdVx = np.delete(rnemdVx, 0)
108 rnemdVy = np.delete(rnemdVy, 0)
109 rnemdZPos = np.delete(rnemdZPos, 0)
110 for i in range(0, numDel+1):
111 rnemdVx = np.delete(rnemdVx, len(rnemdVx)-1)
112 rnemdVy = np.delete(rnemdVy, len(rnemdVy)-1)
113 rnemdZPos = np.delete(rnemdZPos, len(rnemdZPos)-1)
114
115
116 #curve_fit(function to fit to, x-data, y-data, [initial guess of param 1, initial guess of param 2,...]
117 if (jzType == 0):
118 popt, pcov = curve_fit(func, rnemdZPos, rnemdVx, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
119 elif (jzType == 1):
120 popt, pcov = curve_fit(func, rnemdZPos, rnemdVy, [Vs, Z1, Z2, m, w, Vl], ftol=0.01)
121
122
123 #write out the raw data and the fit for visual inspection of quality of fit
124 print("Writing out fit file")
125 fitFileName = outputFileName + "FIT"
126 fitFile = open(fitFileName, "w")
127 fitFile.write("# zPosition <V_{i}> fit value" + "\n")
128 for i in range(0, len(rnemdZPos)):
129 y = func(rnemdZPos[i], popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
130 if (jzType == 0):
131 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVx[i]) + "\t" + str(y) + "\n")
132 elif (jzType == 1):
133 fitFile.write(str(rnemdZPos[i]) + "\t" + str(rnemdVy[i]) + "\t" + str(y) + "\n")
134
135
136 #print "popt = ", popt
137 return (popt, pcov)
138
139
140
141
142
143
144
145def calcShearRate(jzType, rnemdZPos, rnemdV):
146 # jzType = 0 -> flux was in x-dimension
147 # jzType = 1 -> flux was in y-dimension
148
149 velLiq = (rnemdV[0][jzType] + rnemdV[len(rnemdV)-1][jzType])/2.0
150
151 if (len(rnemdV) % 2) == 0 :
152 velSol = ( rnemdV[(len(rnemdV) / 2)-1][jzType] + rnemdV[(len(rnemdV)/ 2)][jzType]) / 2.0
153 elif (len(rnemdV) % 2) != 0:
154 velSol = ( rnemdV[(len(rnemdV) + 1)/2][jzType] + rnemdV[((len(rnemdV)+1)/2)+1][jzType]) / 2.0
155
156 shearRate = velSol - velLiq
157
158 return (shearRate)
159
160
161
162
163
164
165def lambdaCalc(rnemdJzP, jzType, popt):
166 print("Performing lambda calculation")
167 #popt = Vs, Z1, Z2, ml, w, Vl
168
169 zLiq = (popt[0] - popt[5]) / popt[3] + popt[1]
170 zSol = popt[1] + 0.5*popt[4]
171 delta = zLiq - zSol
172
173 if (jzType == 0):
174 lambda1 = rnemdJzP[0] / (delta * popt[3])
175 elif (jzType == 1):
176 lambda1 = rnemdJzP[1] / (delta * popt[3])
177
178
179 return (zLiq, zSol, delta, lambda1)
180
181
182
183
184
185
186
187
188
189
190
191def writeOutputFile(outputFileName, rnemdFileName, popt):
192 outFile = open(outputFileName, "w")
193
194 outFile.write("##################################################### \n")
195 outFile.write("## This output file was generated by slipLength on " + str(datetime.datetime.now()) + " ## \n#\n")
196 outFile.write("# The velocity profile found in " + str(rnemdFileName) + " was fit by \n")
197 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")
198 outFile.write("# where (Vs, Vl, Z1, Z2, ml, w) are fit parameters, and here y = V(i) and x = z \n")
199 outFile.write("# Vs = velocity of the solid\n# Vl = velocity of the liquid at the interface \n")
200 outFile.write("# Z1 = z-position of the lower interface \n# Z2 = z-position of the upper interface \n")
201 outFile.write("# ml = slope of the velocity profile in the lower liquid region of the box \n")
202
203 outFile.write("# w = width of the interface \n#\n")
204 outFile.write("# Obtained optomized parameters \n")
205 outFile.write("# \t Vs = " + str(popt[0]) + "\n")
206 outFile.write("# \t Vl = " + str(popt[5]) + "\n")
207 outFile.write("# \t Z1 = " + str(popt[1]) + "\n")
208 outFile.write("# \t Z2 = " + str(popt[2]) + "\n")
209 outFile.write("# \t ml = " + str(popt[3]) + "\n")
210 outFile.write("# \t w = " + str(popt[4]) + "\n#\n")
211
212 if (float(popt[0]) < float(popt[5])):
213 print("Bad fit detected, Vs < Vl")
214 outFile.write("# WARNING Vs < Vl \n#\n")
215
216
217
218
219
220def writeOutputLambda(outputFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate):
221 outFile = open(outputFileName, "a")
222 outFile.write("# The interfacial slip length, \lambda, was calculated by \n")
223 outFile.write("# \tlambda = Jz(p) / (dV/dz)_{liquid} / delta \n")
224 outFile.write("# where Jz(p) is the imposed momentum flux, (dV/dz)_{liquid} is the slope of the velocity profile in the liquid \n")
225 outFile.write("# portion of the simulation box, and delta is the slip length determined by projecting the velocity profile of \n")
226 outFile.write("# liquid into the solid. \n#\n")
227 outFile.write("# zLiq = " + str(zLiq) + "\n")
228 outFile.write("# zSol = " + str(zSol) + "\n")
229 outFile.write("# \delta (angstroms) = " + str(delta) + "\n")
230 if (jzType == 0):
231 outFile.write("# Jz(px) = " + str((float(rnemdJzP[0]))) + "\n#\n")
232 elif (jzType == 1):
233 outFile.write("# Jz(py) = " + str((float(rnemdJzP[1]))) + "\n#\n")
234 outFile.write("# shear rate (angstroms / fs) = " + str((float(shearRate))) + "\n" )
235 outFile.write("# \lambda (amu/fs/angstrom/angstrom) = " + str((float(lambda1))) + "\n")
236
237 if (lambda1 < 0.0):
238 outFile.write("# WARNING NEGATIVE LAMBDA VALUE FOUND \n#\n#\n")
239
240
241
242
243
244
245
246
247
248
249
250
251def main(argv):
252 parser = argparse.ArgumentParser(
253 description='OpenMD solid/liquid slip length and slip-boundary friction coefficient calculator for orthorhombic systems.',
254 #formatter_class=RawDescriptionHelpFormatter,
255 epilog="Example: slipLength -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")
256 parser.add_argument("-i", "--rnemd-file=", action="store", dest="rnemdFileName", help="use specified input (.rnemd) file")
257 parser.add_argument("-o", "--vfit-file=", action="store", dest="vfitFileName", help="use specified output (.vfit) file")
258 parser.add_argument("-z1", "--lowerGibbsZ=", action="store", type=float, dest="z1", help="the location of the lower Gibbs dividing surface")
259 parser.add_argument("-z2", "--upperGibbsZ=", action="store", type=float, dest="z2", help="the location of the upper Gibbs dividing surface")
260 parser.add_argument("-l", "--lowerZVal=", action="store", nargs='?', type=float, dest="l", help="the initial estimate of the lower interface location (default=z1)")
261 parser.add_argument("-u", "--upperZVal=", action="store", nargs='?', type=float, dest="u", help="the initial estimate of the upper interface location (default=z2)")
262 parser.add_argument("-w", "--intWidth=", action="store", type=float, dest="w", help="the width of the interface")
263 parser.add_argument("-Vs", "--solidVel=", action="store", type=float, dest="Vs", help="the initial estimate of the velocity of the solid")
264 parser.add_argument("-Vl", "--liquidVel=", action="store", type=float, dest="Vl", help="the initial estimate of the velocity of the liquid")
265 parser.add_argument("-m", "--liquidSlope=", action="store", type=float, dest="m", help="the initial estimate of the slope in the liquid")
266 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)")
267
268
269
270
271 if len(sys.argv) == 1:
272 parser.print_help()
273 sys.exit(2)
274 args = parser.parse_args()
275
276 if (not args.rnemdFileName):
277 parser.print_help()
278 parser.error("No input-file was specified")
279
280 if (not args.vfitFileName):
281 parser.print_help()
282 parser.error("No output-file was specified")
283
284 if (not args.z1):
285 parser.print_help()
286 parser.error("No lower Gibbs dividing surface specified")
287
288 if (not args.z2):
289 parser.print_help()
290 parser.error("No upper Gibbs dividing surface specified")
291
292 if (not args.l):
293 args.l = args.z1
294
295 if (not args.u):
296 args.u = args.z2
297
298 if (not args.w):
299 parser.print_help()
300 parser.error("No width of interface specified")
301
302 if (not args.Vs):
303 parser.print_help()
304 parser.error("No initial estimate of the solid velocity specified")
305
306 if (not args.Vl):
307 parser.print_help()
308 parser.error("No initial estimate of the liquid velocity specified")
309
310 if (not args.m):
311 parser.print_help()
312 parser.error("No initial estimate of the slope in the liquid velocity profile specified")
313
314
315
316
317
318
319 #Call functions here, pass appropriate variables.
320
321
322 (rnemdJzP, jzType, rnemdZPos, rnemdV) = rnemdExtractor(args.rnemdFileName)
323 (popt, pcov) = velFitter(rnemdJzP, jzType, rnemdZPos, rnemdV, args.Vs, args.l, args.u, args.m, args.w, args.Vl, args.toDel, args.vfitFileName)
324
325
326
327 (shearRate) = calcShearRate(jzType, rnemdZPos, rnemdV)
328 (zLiq, zSol, delta, lambda1) = lambdaCalc(rnemdJzP, jzType, popt)
329
330
331
332 writeOutputFile(args.vfitFileName, args.rnemdFileName, popt)
333 writeOutputLambda(args.vfitFileName, rnemdJzP, jzType, zLiq, zSol, delta, lambda1, shearRate)
334
335
336 print("Summary of calculation")
337 if (jzType == 0):
338 print("Jz(p) = ", float(rnemdJzP[0]))
339 elif (jzType == 1):
340 print("Jz(p) = ", float(rnemdJzP[1]))
341
342 print("Shear rate (angstroms / fs) = ", float(shearRate))
343 print("Slip length \delta (angstroms) = ", float(delta))
344 print("Friction coefficient \lambda (amu / fs / angstroms^2) = ", float(lambda1))
345
346
347if __name__ == "__main__":
348 main(sys.argv[1:])