#making groups variable r equal 1 variable a equal 0.0 shell mkdir dir$r/dump shell mkdir dir$r/restart shell mkdir dir$r/logfile shell mkdir dir$r/output shell mkdir dir$r/output_step variable conc equal 1 variable sigma equal 19.79 # Gaussian parameter in nm?¹ variable q0 equal 0.0 variable total_charge_per_atom equal v_q0*erfc(v_sigma*sqrt(pi)*0.5) # Define the Gaussian charge distribution using fix addforce ###################################### units real atom_style full dimension 3 newton on processors * * 2 boundary p p f read_data input_data/embf240.dat #read_restart input_data/x.r_1.39000 ################################################ ############force field######################## ############################################ pair_style lj/cut/coul/long 14 bond_style harmonic angle_style harmonic dihedral_style opls kspace_style pppm/electrode 1.0e-6 kspace_modify slab 3 fix zwalls all wall/reflect zlo EDGE zhi EDGE units box pair_coeff 1 1 1.12 5.06 # PF6(A) pair_coeff 2 2 0.0474 3.211 # A1 pair_coeff 3 3 0.0474 3.211 # A2 pair_coeff 4 4 0.1 3.40 # C pair_coeff 5 5 0.612 4.38 # C1 pair_coeff 6 6 0.86 3.41 # C2 pair_coeff 7 7 0.29 4.38 # C3 pair_coeff 8 8 0.38 3.60 # Me pair_coeff 9 9 0.1 3.30 # N #pair done pair_modify shift yes mix arithmetic bond_coeff 1 385.0 1.458 bond_coeff 2 385.0 1.458 bond_coeff 3 385.0 1.458 bond_coeff 4 650.0 1.157 #C-N #bonds done angle_coeff 1 70 153.94 # c2-c1-c3 angle_coeff 2 70 178.6 # ME-C-N #angle done #dihedral_modify mix arithmetic mass 1 144.96 # A mass 2 12.0107 # A1 mass 3 12.0107 # A2 mass 4 12.0107 # C mass 5 67.07 # C1 mass 6 15.04 # C2 mass 7 29.07 # C3 mass 8 15.04 # Me mass 9 14.0067 # N #systems group emim type 5 6 7 group pf6 type 1 #group bf4 type 1 group elec1 type 2 group elec2 type 3 group sol type 1 4 5 6 7 8 9 group syst type 1 2 3 4 5 6 7 8 9 group ele type 2 3 group elel type 2 group eler type 3 group elelr molecule 44444 group elell molecule 44445 group elerr molecule 44441 group elerl molecule 44442 group acn type 4 8 9 #charges set type 1 charge -0.78 # BF4 or PF6 set type 2 charge 0 # A1 set type 3 charge 0 # A2 set type 4 charge 0.129 # C set type 5 charge 0.3591 # C1 set type 6 charge 0.1888 # C2 set type 7 charge 0.2321 # C3 set type 8 charge 0.269 # Me set type 9 charge -0.398 # N #charges done #modification #fix charge all addforce ${total_charge_per_atom}*erfc(v_sigma*sqrt(pi)*(x-x0))*exp(-v_sigma^2*(x-x0)^2) run_style verlet neighbor 2.0 bin neigh_modify exclude group ele ele check no neigh_modify delay 0 every 1 check yes log dir$r/logfile/log_${conc} timestep 0.05 thermo 100 compute t_sol sol temp fix 5 ele setforce 0.0 0.0 0.0 fix 6 sol shake 0.0001 20 0 b 1 minimize 1.0e-6 1.0e-8 100 10000 velocity sol create 298.0 86289 timestep 1 dump 1 syst custom 100 dir$r/dump/dump-${conc}.lammpstrj id mol type x y z q dump_modify 1 sort id dump 2 syst custom 5000 dir$r/dump/dumpall-${conc}.lammpstrj id mol type x y z q dump_modify 2 sort id restart 100000 dir$r/output/x.r_${conc} restart 1000 dir$r/output_step/x.r_${conc} shell echo "Hello, this is a message in the LAMMPS log." dump 4 syst custom 100 dir$r/dump/dump1.lammpstrj id type x y z q dump_modify 4 sort id compute emsd emim msd com yes variable twopoint1 equal c_emsd[4]/4/(step*dt+1.0e-6) fix 9 all vector 10 c_emsd[4] variable fitslope1 equal c_emsd[4]/4/(dt) compute pmsd pf6 msd com yes variable twopoint2 equal c_pmsd[4]/4/(step*dt+1.0e-6) fix 10 all vector 10 c_pmsd[4] variable fitslope2 equal c_pmsd[4]/4/(dt) region region1 block 0.0375 50.0375 0.0035 50.0035 -100.00 100.00 units box compute emim_num_density emim chunk/atom bin/1d z -100.00 2 region region1 units box compute pf_num_density pf6 chunk/atom bin/1d z -100.00 2 region region1 units box compute acn_num_density acn chunk/atom bin/1d z -100.00 2 region region1 units box fix num_em emim ave/chunk 100 1 100 emim_num_density density/number ave one file num_emim_density.txt fix num_pf pf6 ave/chunk 100 1 100 pf_num_density density/number ave one file num_pf6_density.txt fix num_acn acn ave/chunk 100 1 100 acn_num_density density/number ave one file num_acn_density.txt compute overall_density sol chunk/atom bin/1d z -102.00 1 region region1 units box fix dens_overall sol ave/chunk 1000 1 1000 overall_density density/mass ave one file densities.txt #compute emim_COM emim com yes #fix 1 emim ave/time 100 1 100 c_emim_COM[*] file emimcom.txt mode vector timestep 0.01 fix 3 sol nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_sol run 50000 unfix 3 write_restart dir$r/finalpart1.rest write_data dir$r/finalpart1.dat timestep 0.1 fix 3 sol nvt temp 298.0 298.0 10000.0 fix_modify 3 temp t_sol run 50000 unfix 3 write_restart dir$r/finalpart1.rest write_data dir$r/finalpart1.dat timestep 1 fix 3 sol nvt temp 298.0 298.0 10000.0 fix_modify 3 temp t_sol run 100000 unfix 3 timestep 2 fix 3 sol nvt temp 298.0 298.0 10000.0 fix_modify 3 temp t_sol run 200000 unfix 3 timestep 2 variable w equal ${a}/2 fix e elec1 electrode/conp -${w} 1.979 couple elec2 ${w} symm on etypes on write_inv /home/kuntekaushik/33kcg2.csv fix 4 sol nvt temp 298.0 298.0 10000.0 fix_modify 4 temp t_sol variable q atom q compute ql elel reduce sum v_q compute qr eler reduce sum v_q compute qll elell reduce sum v_q compute qlr elelr reduce sum v_q compute qrl elerl reduce sum v_q compute qrr elerr reduce sum v_q thermo_style custom step c_ql c_qr c_qll c_qlr c_qrl c_qrr pe v_twopoint1 v_fitslope1 v_twopoint2 v_fitslope2 thermo_modify flush yes norm no run 1000000