#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/embf15M300.dat #read_restart input_data/x.r_1.624000 ################################################ ############force field######################## ############################################ pair_style lj/cut/coul/long 12 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 0.775 4.51 # BF4(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 350.0 1.3 bond_coeff 3 350.0 1.3 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 86.81 # 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 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 elell molecule 44444 group elelr molecule 44445 group elerl molecule 66665 group elerr molecule 66666 group elel molecule 44444 44445 group eler molecule 66666 66665 group acn type 4 8 9 #charges set type 1 charge -0.78 # BF4 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 1000 compute t_all all temp fix 5 ele setforce 0.0 0.0 0.0 fix 6 sol shake 0.0001 20 0 b 4 minimize 1.0e-6 1.0e-8 150 10000 velocity sol create 298.0 86289 dump 1 syst custom 500 dir$r/dump/dump-${conc}.lammpstrj id mol type x y z q dump_modify 1 sort id dump 2 syst custom 20000 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." ############################################################### #compute dens sol chunk/atom bin/1d z lower 0 upper 100 0.01 units reduced #fix 7 sol ave/chunk 10 1 100 dens density/mass ave window 10 file densities.txt #region my_region block -0.0195 49.9805 -0.028499 49.971501 0 100 #group my_group region my_region #variable atoms_in_region equal count(my_group) # Print the count to the log file #print "Number of atoms in the region: ${atoms_in_region}" #compute my_temp my_group temp #thermo_modify temp my_temp #thermo_style custom step c_temp ############################################################ #thermo_style custom step c_my_temp dump 4 syst custom 1000 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 bmsd bf4 msd com yes variable twopoint2 equal c_bmsd[4]/4/(step*dt+1.0e-6) fix 10 all vector 10 c_bmsd[4] variable fitslope2 equal c_bmsd[4]/4/(dt) region region1 block -0.0025 50.0025 0.002001 50.002001 -102 102 units box compute emim_num_density emim chunk/atom bin/1d z -102 2 region region1 units box compute bf_num_density bf4 chunk/atom bin/1d z -102 2 region region1 units box compute acn_num_density acn chunk/atom bin/1d z -102 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_bf bf4 ave/chunk 100 1 100 bf_num_density density/number ave one file num_bf4_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 -51.999001 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 all nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_all run 10000 unfix 3 timestep 0.1 fix 3 all nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_all run 10000 unfix 3 timestep 1 fix 3 all nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_all run 20000 unfix 3 timestep 2 fix 3 all nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_all 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/embf1533kcg.csv fix 4 all nvt temp 298.0 298.0 100 fix_modify 4 temp t_all 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