# Input file for RBCh flow test label rerun variable latden index 3 0 variable datafile string rbc_D50 # Sets the data file name variable P_damp equal 1 # Sets the Pressure damping variable used in the following " fix npt" commands variable P equal 0.04 #0.05 # Sets the desired pressure of the system used in the following " fix npt" commands variable w index 2.3 2.7 0 #2.7 0 variable j index 3.3 3.5 0 variable scale1 equal ramp(2.7,$w) # Sets the scale performed in the final fix on water scaling down sigma from 2.7- 2.5 variable scale2 equal ramp(2.7,$j) variable ew equal 0.2 # Sets the energy well ( epsilon) variable used in the pair coefficient specifications variable ini_T equal 0.02 #0.02 # Sets the initial temperature used in the fix commands below variable T equal 0.15 #0.23 # Sets the final temperature used in the fix commands below variable LD equal 1.0 # Sets the temperature damping parameter used in fix commands below variable fid string ${datafile} variable logfid string log.${fid} variable beta equal 0.0194*0.0 variable q index 537 414 225 53 31 69 414 547 97 77 13 62 0 variable f index 4 0 variable a index 75 60 90 105 0 variable flow atom sqrt(40-abs(c_z))*0.001*$f log RBCh_4FS_a$a_f$f_1c$w_2c$j.lammps # ------------------------ INITIALIZATION ---------------------------- units lj dimension 3 boundary p p p atom_style hybrid ellipsoid peri molecular read_data ${datafile}.data extra/atom/types 1 group cell1 mass * 1 #communicate single vel yes comm_modify vel yes # ------------------------ FORCE FIELDS ------------------------------ #pair_style hybrid lj/cut 3.6 fluidmembrane 2.6 # Sets the pair_style potentials used for the simulation pair_style hybrid lj/cut 3.6 fluidmembrane 2.6 dpd 1.1 2.5 $q pair_coeff 1*2 1*2 fluidmembrane 1.0 1.0 2.6 4 3 ${beta} # Sets the pair coefficients between 1 and 2 (bilayer) particles pair_coeff 1*2 3*5 lj/cut ${ew} 1.0 # Sets the pair coefficients between 1,2 (bilayer) and 3,4,5 (network) pair_coeff 1*2 6*8 lj/cut ${ew} 1.0 # Sets the pair coefficients between 1,2 (bilayer) and 6,7 (water) pair_coeff 3*5 3*5 lj/cut ${ew} 1.0 # Sets the pair coefficients between 3,4,5 (network) and network pair_coeff 3*5 6*8 lj/cut ${ew} 1.0 # Sets the pair coefficients between 3,4,5 (network) and 6,7 (water) pair_coeff 6*8 6*8 lj/cut ${ew} 2.7 # Sets the pair coefficients between 6,7 (water) and 6,7 (water) pair_coeff * 9 dpd 5 1 pair_coeff 9 9 lj/cut 0 0 neighbor 1.0 bin neigh_modify check yes every 1 delay 0 page 1000000 one 100000 bond_style harmonic1 bond_coeff 1 50 1.5 # spring constant for the spectrin network bonds bond_coeff 2*3 20 2.0 # spring constant for the transmembrane protein to spectrin network # ---------------------- SPATIAL DEFINITION -------------------------- #change_box all y final -100 100 z final -40 40 x final -75 75 units box change_box all y final -100 100 z final -40 40 x final -100 100 units box read_data ${datafile}.data add append shift -62.7 0 0 group cell2 read_data ${datafile}.data add append shift 0 62.7 0 group cell2 read_data ${datafile}.data add append shift -62.7 62.7 0 group cell1 change_box all x final -100 150 z final -40 40 y final -40 140 units box set type 1*2 shape 1 1 1 # membrane particle are rigid body, it needs shape parameter to be set group bilayer type 1 2 # Creates a group called bilayer which consists of type 1 and 2 particles, #(Type 1 Bilayer, Type 2 Transmembrane proteins) group network type 3 4 5 # Creates a group called network which consists of type 3 4 5 particles group anchor_bonds type 2 3 5 group water_in1 type 6 # Creates a group called water_in consisting of type 6 particles group water_in2 intersect cell2 water_in1 set group water_in2 type 8 group rbc type 1 2 3 4 5 region floor block INF INF INF INF -40 -39 units box region slant1 prism 35 75 -40 120 -20 40 $a 0 0 units box side in region 1cy1a cylinder x -20 0 4 35 150 side in region 1cy1b cylinder x -20 0 3 35 150 side out region 1cy1 intersect 3 1cy1a 1cy1b slant1 region 1cy2a cylinder x -20 10 4 35 150 side in region 1cy2b cylinder x -20 10 3 35 150 side out region 1cy2 intersect 3 1cy2a 1cy2b slant1 region 1cy3a cylinder x -20 20 4 35 150 side in region 1cy3b cylinder x -20 20 3 35 150 side out region 1cy3 intersect 3 1cy3a 1cy3b slant1 region 1cy4a cylinder x -20 30 4 35 150 side in region 1cy4b cylinder x -20 30 3 35 150 side out region 1cy4 intersect 3 1cy4a 1cy4b slant1 region 1cy union 4 1cy1 1cy2 1cy3 1cy4 region 1cyb intersect 4 1cy1b 1cy2b 1cy3b 1cy4b region 2cy1a cylinder x 10 0 4 35 150 side in region 2cy1b cylinder x 10 0 3 35 150 side out region 2cy1 intersect 3 2cy1a 2cy1b slant1 region 2cy2a cylinder x 10 10 4 35 150 side in region 2cy2b cylinder x 10 10 3 35 150 side out region 2cy2 intersect 3 2cy2a 2cy2b slant1 region 2cy3a cylinder x 10 20 4 35 150 side in region 2cy3b cylinder x 10 20 3 35 150 side out region 2cy3 intersect 3 2cy3a 2cy3b slant1 region 2cy4a cylinder x 10 30 4 35 150 side in region 2cy4b cylinder x 10 30 3 35 150 side out region 2cy4 intersect 3 2cy4a 2cy4b slant1 region 2cy union 4 2cy1 2cy2 2cy3 2cy4 region 2cyb intersect 4 2cy1b 2cy2b 2cy3b 2cy4b region 3cy1a cylinder x 40 0 4 35 150 side in region 3cy1b cylinder x 40 0 3 35 150 side out region 3cy1 intersect 3 3cy1a 3cy1b slant1 region 3cy2a cylinder x 40 10 4 35 150 side in region 3cy2b cylinder x 40 10 3 35 150 side out region 3cy2 intersect 3 3cy2a 3cy2b slant1 region 3cy3a cylinder x 40 20 4 35 150 side in region 3cy3b cylinder x 40 20 3 35 150 side out region 3cy3 intersect 3 3cy3a 3cy3b slant1 region 3cy4a cylinder x 40 30 4 35 150 side in region 3cy4b cylinder x 40 30 3 35 150 side out region 3cy4 intersect 3 3cy4a 3cy4b slant1 region 3cy union 4 3cy1 3cy2 3cy3 3cy4 region 3cyb intersect 4 3cy1b 3cy2b 3cy3b 3cy4b region 4cy1a cylinder x 70 0 4 35 150 side in region 4cy1b cylinder x 70 0 3 35 150 side out region 4cy1 intersect 3 4cy1a 4cy1b slant1 region 4cy2a cylinder x 70 10 4 35 150 side in region 4cy2b cylinder x 70 10 3 35 150 side out region 4cy2 intersect 3 4cy2a 4cy2b slant1 region 4cy3a cylinder x 70 20 4 35 150 side in region 4cy3b cylinder x 70 20 3 35 150 side out region 4cy3 intersect 3 4cy3a 4cy3b slant1 region 4cy4a cylinder x 70 30 4 35 150 side in region 4cy4b cylinder x 70 30 3 35 150 side out region 4cy4 intersect 3 4cy4a 4cy4b slant1 region 4cy union 4 4cy1 4cy2 4cy3 4cy4 region 4cyb intersect 4 4cy1b 4cy2b 4cy3b 4cy4b region 5cy1a cylinder x 100 0 4 35 150 side in region 5cy1b cylinder x 100 0 3 35 150 side out region 5cy1 intersect 3 5cy1a 5cy1b slant1 region 5cy2a cylinder x 100 10 4 35 150 side in region 5cy2b cylinder x 100 10 3 35 150 side out region 5cy2 intersect 3 5cy2a 5cy2b slant1 region 5cy3a cylinder x 100 20 4 35 150 side in region 5cy3b cylinder x 100 20 3 35 150 side out region 5cy3 intersect 3 5cy3a 5cy3b slant1 region 5cy4a cylinder x 100 30 4 35 150 side in region 5cy4b cylinder x 100 30 3 35 150 side out region 5cy4 intersect 3 5cy4a 5cy4b slant1 region 5cy union 4 5cy1 5cy2 5cy3 5cy4 region 5cyb intersect 4 5cy1b 5cy2b 5cy3b 5cy4b region cyb intersect 5 1cyb 2cyb 3cyb 4cyb 5cyb region cy union 5 1cy 2cy 3cy 4cy 5cy region slant11 intersect 2 slant1 cyb region slant2 prism 36 74 -39 119 -19 39 $a 0 0 units box side out region slant3 intersect 2 slant11 slant2 region slant union 2 slant3 cy lattice fcc ${latden} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 create_atoms 9 region slant create_atoms 9 region floor group slant type 9 #region ww block -90 50 40 90 -20 20 units box #lattice sc 0.06 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 #create_atoms 7 region ww group water_out type 7 # Creates a group called water_out consisting of type 7 particles group water type 6 7 8 # Creates a group called water consisting of both type 6 and 7 particles # ------------------------- SETTINGS --------------------------------- compute peatom all pe/atom compute peall all pe compute keatom all ke/atom compute keall all ke compute temp all temp compute astr rbc stress/atom temp compute x all property/atom x compute xm all reduce ave c_x compute y all property/atom y compute ym all reduce ave c_y compute z all property/atom z compute zm all reduce ave c_z compute xl all reduce min c_x compute yl all reduce min c_y compute zl all reduce min c_z compute xh all reduce max c_x compute yh all reduce max c_y compute zh all reduce max c_z compute vx all property/atom vx compute vy all property/atom vy compute vz all property/atom vz compute fx all property/atom fx compute fy all property/atom fy compute fz all property/atom fz variable nall equal count(all) variable peall equal c_peall variable keall equal c_keall variable peatom atom c_peatom variable keatom atom c_keatom variable astress atom c_astress variable lstr atom sqrt(((c_astr[1]-c_astr[2])^2+(c_astr[2]-c_astr[3])^2+(c_astr[1]-c_astr[3])^2+6*(c_astr[4]^2+c_astr[5]^2+c_astr[6]^2))/2) compute lstress rbc reduce sum v_lstr variable lstress equal c_lstress ###################################### #timestep 0.01 #velocity all create 1.1 12$q mom yes rot yes velocity bilayer create ${T} 872$q loop geom # Creates Initial Velocity for bilayer particles ( Type 1 and 2) velocity network create ${ini_T} 887$q loop geom # Creates Initial Velocity for network particles ( Type 3 and 4 and 5) velocity water create ${T} 28$q loop geom # Creates Initial Velocity for water particles ( Type 6 and 7 ) thermo_style custom step temp press ebond v_nall # Printing thermodynamic data to the screen and log file thermo 200 # Sets the number of timesteps to run before printing out thermodynamic data thermo_modify lost warn timestep 0.005 dump 1 all cfg 10000 RBCh_4FS_a$a_f$f_1c$w_2c$j_*.cfg mass type xs ys zs type c_astr[1] c_astr[2] c_astr[3] & c_astr[4] c_astr[5] c_astr[6] v_lstr c_peatom c_keatom c_vx c_vy c_vz c_fx c_fy c_fz dump_modify 1 element P O C N F H H H B dump_modify 1 pad 6 run 0 fix out all ave/time 1 100 100 v_keall v_peall v_lstress file RBCh_4FS_a$a_f$f_1c$w_2c$j.txt ################################ equilibrium ########################################################## fix 1 water npt temp ${T} ${T} ${LD} iso ${P} ${P} ${LD} # Applies an NPT ensemble fix to the water group minimize 1.0e-5 1.0e-5 5000 50000 #reset_timestep 0 run 10000 fix 2 network npt temp ${ini_T} ${T} ${LD} iso ${P} ${P} ${P_damp} # Applies an NPT ensemble fix to network group #fix 2 network nve #fix 2 network nvt temp ${ini_T} ${T} ${LD} run 10000 #fix 3 bilayer nvt/asphere temp ${T} ${T} ${LD} # Applies an NVT to bilayer using only temperature parameters fix 3 bilayer nve/asphere run 10000 dump_modify 1 every 50000 ### change RBC volume fix 11 water adapt 1 pair lj/cut sigma 6 6 v_scale1 # Applies a final "adapt" fix to the water group and this will change the sigma value of LJ potential fix 22 water adapt 1 pair lj/cut sigma 8 8 v_scale2 run 50000 unfix 11 unfix 22 run 10000 dump_modify 1 every 20000 #fix push bilayer addforce v_flow 0.0 0.0 fix push water_out addforce v_flow 0.0 0.0 run 400000 ###################################### # SIMULATION DONE clear next q print "Run Done" next w if $w!=0 then "jump RBCh_4FS.in rerun" variable w delete next j if $j!=0 then "jump RBCh_4FS.in rerun" variable j delete next f if $f!=0 then "jump RBCh_4FS.in rerun" variable f delete next a if $a!=0 then "jump RBCh_4FS.in rerun" variable a delete print "Finally Finished" #######################################################################