#*************************************************************************************************** # Liquid argon around CNT #*************************************************************************************************** # Purpose: 1. Setup argon around cnt in equlibrium. # 2. Then impose a uniform flow to study drag oon cnt ( 2 references below ) #-------------------------------------------------------------------------------------------------- # Reference: # # 1. "Drag on a nanotube in uniform liquid argon flow" # W.Tang and S.G.Advani, # Dept. of Mechanical Engineering, Univ. of Delaware. # J. Chemical Physics, vol. 125, 2006 # # 2. "Molecular dynamics simulations of liquid flow in and around carbon nanotubes" # W.D.Nicholls, M.K.Borg and J.M.Reese # Dept. of Mechanical Engineering, Univ. of Strathclyde, Glasgow, UK. # ASME 2010, FEDSM-ICNMM 2010-30360 # (This paper uses paper 1 for reference also ) # #*************************************************************************************************** # Set-up geometry and atoms #*************************************************************************************************** # Basic simulation parameters from reference 1 # # 1. Liquid argon: T = 95 Kelvin ; Density = 1342 kg/m^3 # # 2. CNT: (12,0)-zigzag ; Diameter 9.422 angstroms ; Length = 21.3 angs. # # 3. Domain size: x = 305.6 angs ( flow direction ; flow normal to cnt axis ) # y = 254.7 angs ( cross-flow direction ) # z = 21.3 angs ( channel width -- also length of CNT and cnt axis ) #--------------------------------------------------------------------------------------------------- # Initialization # -------------- units real # "lj, real, SI, ...." dimension 3 # 3 Dimensional flow domain boundary p p p # periodic bc in x,y,z directions atom_style atomic # atom type is "atomic, molecular, full...." # Read CNT data from VMD output #------------------------------ read_data cnt_atomic.data #( Modified inside data file: # 1. atom type 2. domain size 3. add argon atom mass ) # Geometry #--------- # Length dimensions are in angstroms variable xmax_box equal 304.60 variable ymax_box equal 254.70 variable zmax_box equal 24.106 # ( longer than in reference paper ) # Create a rectangular channel region # ID shape x_min x_max y_min y_max z_min z_max key arg region channel block 0 ${xmax_box} 0 ${ymax_box} 0 ${zmax_box} units box # Simulation domain is this channel region ( size defined inside read_data file ) # Create a cylindrical hole larger than cnt # and fill the region outside with argon atoms # so that there are no argon inside cnt or overlaps cnt. # --------------------------------------------------------------------------- # Cylinder geometrical parameters # ------------------------------- variable cyl_x_cen equal 152.30 # cylinder: x coordinate of center variable cyl_y_cen equal 127.35 # cylinder: y coordinate of center variable cyl_rad equal 10.00 # cylinder: radius variable cyl_zlo equal 0.00 # cylinder: length - lower boundary variable cyl_zhi equal v_zmax_box # cylinder: length - upper boundary # ID shape axis x-center y-center radius z-lo z-hi key arg key arg region cyl_out cylinder z ${cyl_x_cen} ${cyl_y_cen} ${cyl_rad} ${cyl_zlo} ${cyl_zhi} side out units box # create a region with a hole at the center of the channel # -------------------------------------------------------- # ID style # of region region key arg # regions ID-1 ID-2 region channel_hole intersect 2 channel cyl_out units box # Fill the region with argon atoms #--------------------------------- # Number of argon atoms for a chosen density variable N_argon equal 37835 # 33670 # Atom-type style #atom rnd-seed region-ID key arg create_atoms 2 random ${N_argon} 3495768 channel_hole units box #*************************************************************************************************** # GROUPING ATOMS # --------------- # group-ID style style arguments group cnt type 1 # create cnt : group all carbon atoms into "CNT" # Move cnt from (x,y)=(0,0) to center of domain. (VMD creates cnt at (0,0) ) # group-ID style dx dy dz key arg displace_atoms cnt move ${cyl_x_cen} ${cyl_y_cen} 0.0 units box # group-ID style style arguments group flow type 2 # create flow: group all argon atoms for flow #*************************************************************************************************** # Channel volume ( ignoring cnt volume ) # -------------------------------------- variable x_box_len equal v_xmax_box variable y_box_len equal v_ymax_box variable z_box_len equal v_zmax_box variable vol_channel equal " v_x_box_len*v_y_box_len*v_z_box_len " #*************************************************************************************************** # Atom properties #*************************************************************************************************** # Mass of atoms #-------------- # Atom type mass value mass 1 12.011 # molar mass of Carbon ; grams/mol mass 2 39.948 # molar mass of Argon ; grams/mol variable mass_gm_mol_carbon equal 12.011 # carbon mass: gm/mol variable mass_gm_mol_argon equal 39.948 # argon mass: gm/mol # Atom potential #--------------- #[ Carbon in cnt are fixed, and "atomic_style" ] variable ljcutoff equal 7.65 # cut-off: ; Angstroms variable eps_aa equal 0.2375 # epsilon: argon--argon ; Kcal/mol variable sig_aa equal 3.405 # sigma: argon--argon ; Angstroms variable eps_ac equal 0.2816 # epsilon: argon--carbon ; Kcal/mol variable sig_ac equal 3.573 # sigma: argon--carbon ; Angstroms # potential-type potential-type arguments pair_style lj/cut ${ljcutoff} # LJ potential ; cut-off size # Atom-type potential parameters # Interaction pair_coeff 1 1 ${eps_ac} ${sig_ac} # LJ potential: carbon--carbon (dummy) pair_coeff 1 2 ${eps_ac} ${sig_ac} # LJ potential: argon--carbon pair_coeff 2 2 ${eps_aa} ${sig_aa} # LJ potential: argon--argon # eps: kcal/mol sigma: Angstroms #*************************************************************************************************** # Published Argon potentials and characteristic scales : # Variable symbol value units # ---------------------------------------------------- # 1. sigma s 0.3405 nm (nano-meters) ; 1 nm = 10 Angstrom # 3.405 Angstroms # 2. epsilon e 119.8*K_b Joules ; K_b = Boltzmann const. in J/K # 1.65*10(-21) Joules # 0.2375 kcal/mol ; 1 cal = 4.184 J . # 3. molar mass 39.948 gram/mole # 4. mass (atom) 6.63*10(-26) kg # 6.69*10(-26) kg # 5. Time 2.17*10(-12) seconds ; 1 sec = 10(15) femto-sec # 2170 femtoseconds # 6. Velocity 157 meters/second ; 1 m = 10(10) Angstrom # 0.00157 angstrom/femtosec # 7. Temperature e/k_b 119.8 Kelvin #*************************************************************************************************** # Density of argon in channel # Conversion factors variable avogadro equal 6.022*(10^(23)) variable conv_gm_kg equal 0.001 variable conv_angs_m equal 10^(-10) variable conv_atm_pascal equal 101325 variable mass_atom_gm_argon equal ${mass_gm_mol_argon}/${avogadro} variable mass_atom_kg_argon equal ${mass_atom_gm_argon}*${conv_gm_kg} # Argon mass variable argon_ch_num_atom equal count(flow) variable argon_mass equal ${mass_atom_kg_argon}*${argon_ch_num_atom} # Argon density variable vol_channel_m equal ${vol_channel}*((${conv_angs_m})^3) variable argon_density equal ${argon_mass}/${vol_channel_m} #--------------------------------------------------------------------------------------------------- # 1. Argon density print " Argon density Kg/m3 = ${argon_density} " # 2. Count carbon atoms variable cnt_num equal count(cnt) print " CNT: Number of atoms = ${cnt_num} " #quit #*************************************************************************************************** # Velocity of atoms # ----------------- # variable-name style variable-value variable T equal 95.00 # set inital temperature variable rndG equal 3584978 # set initial random seed # group-ID style style arguments key arg key arg velocity flow create $T ${rndG} dist gaussian units box # Gaussian dist. #*************************************************************************************************** # Computing thermodynamic variables # --------------------------------- # compute-ID group-ID compute_style compute therm_channel flow temp # temperature compute pr_channel all pressure therm_channel # pressure # style style-arg: thermodynamic variables from "compute" c_* thermo_style custom c_therm_channel c_pr_channel etotal # store/print (T,P.E) # dummy run to evaluate thermodynamic variables in "compute" commands run 0 # dummy run to generate thermal variables # Assign thermal variables value from dummy run 0 #------------------------------------------------ # variable-name style variable-value variable temper_channel equal c_therm_channel # thermal variables from thermo_style memory variable press_channel equal c_pr_channel # Pressure in kilo-pascals # [ to check with reference paper; lammps pressure is in atmos. for "real" units] variable press_channel_kpa equal (v_press_channel*${conv_atm_pascal})/1000.0 print "Temperature Pressure: Channel: = ${temper_channel} ${press_channel_kpa}" #*************************************************************************************************** # Methods to remove atom overlap from random creation #*************************************************************************************************** #--------------------------------------------------------------------------------------------------- # Method 1 : Use soft potential to push apart overlapping atoms #--------------------------------------------------------------------------------------------------- # Use same parameters as argon-argon potential pair_style soft 3.405 # cut-off == sigma : argon-argon pair_coeff * * 0.2375 # epsilon: argon-argon thermo 100 timestep 4 fix nve_smooth flow nve run 10000 unfix nve_smooth #quit #--------------------------------------------------------------------------------------------------- # Method 2 : Energy minimization to remove atom overlap #--------------------------------------------------------------------------------------------------- # fix cnt atoms during minimization fix cnt_fix cnt setforce 0.0 0.0 0.0 # Energy minimization to relax atoms-overlap due to random atom creation # etol ftol max-iter max-eval minimize 1.0e-4 1.0e-6 1000 10000 # Remove the fix cnt_fix unfix cnt_fix #quit #--------------------------------------------------------------------------------------------------- # Method 3 : Deleting atoms #--------------------------------------------------------------------------------------------------- #variable N_argon_bd equal count(flow) #print " Before delete: = ${N_argon_bd} " #delete_atoms overlap 1.0 flow flow #variable N_argon_ad equal count(flow) #print " After delete: = ${N_argon_ad} " #quit # #--------------------------------------------------------------------------------------------------- # Method 4 : "fix nve/limit" and "fix viscous" #--------------------------------------------------------------------------------------------------- #thermo 10 #timestep 4 #fix vis_smooth flow viscous 1.0 # viscous damping to drain energy #fix nve_smooth flow nve/limit 1.0 # Limit atom displacement over a timestep #run 1000 #unfix vis_smooth #unfix nve_smooth #quit #*************************************************************************************************** # SIMULATION #*************************************************************************************************** #--------------------------------------------------------------------------------------------------- # Reset atom potential after soft potential run. #--------------------------------------------------------------------------------------------------- # # potential-type potential-type arguments pair_style lj/cut ${ljcutoff} # LJ potential ; cut-off size # Atom-type potential parameters # Interaction pair_coeff 1 1 ${eps_ac} ${sig_ac} # LJ potential: carbon--carbon (dummy) pair_coeff 1 2 ${eps_ac} ${sig_ac} # LJ potential: argon--carbon pair_coeff 2 2 ${eps_aa} ${sig_aa} # LJ potential: argon--argon # eps: kcal/mol sigma: Angstroms #--------------------------------------------------------------------------------------------------- # Reset velocity after smoothing inital atom distribution #--------------------------------------------------------------------------------------------------- # group-ID style style arguments velocity flow create $T ${rndG} dist gaussian units box #*************************************************************************************************** # Test run a few steps of NVT simulation #*************************************************************************************************** #--------------------------------------------------------------------------------------------------- # Write atom data out to file to visualize in VMD # ----------------------------------------------- # dump-ID group-ID style Freq. out-file name style arguments # # Time-step dump 1 all custom 100 out_data.custom id type mass x y z vx vy vz #--------------------------------------------------------------------------------------------------- # Simulation # ---------- variable timestep_nvt_damp equal 25 # time steps: Noose-Hoover damp/relax in NVT thermo 10 # output thermal prop every X time steps timestep 4 # time step in femto seconds # fix-ID group-ID fix nvt_min flow nvt temp $T $T ${timestep_nvt_damp} # NVT : constant N, volume, temperature # flow: == time-stepping argon only, keeping cnt fixed. run 1000 unfix nvt_min #quit #--------------------------------------------------------------------------------------------------- #*************************************************************************************************** # END PROGRAM #***************************************************************************************************