#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/embf240f.dat #read_restart dir1/output_step/x.r_1.60000 #read_restart dir$r/finalpart1.rest ################################################ ############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.095 3.5814 # B1 pair_coeff 2 2 0.085 3.55 # C1 pair_coeff 3 3 0.085 3.55 # C2 pair_coeff 4 4 0.085 3.55 # C3 pair_coeff 5 5 0.109 3.55 # C4 pair_coeff 6 6 0.109 3.55 # C5 pair_coeff 7 7 0.109 3.55 # C6 pair_coeff 8 8 0.06 3.1181 # F1 pair_coeff 9 9 0.06 3.1181 # F2 pair_coeff 10 10 0.06 3.1181 # F3 pair_coeff 11 11 0.06 3.1181 # F4 pair_coeff 12 12 0.0474 3.211 # A1 pair_coeff 13 13 0.03 2.42 # H1 pair_coeff 14 14 0.03 2.5 # H10 pair_coeff 15 15 0.03 2.42 # H11 pair_coeff 16 16 0.03 2.42 # H2 pair_coeff 17 17 0.03 2.5 # H3 pair_coeff 18 18 0.03 2.5 # H4 pair_coeff 19 19 0.03 2.5 # H5 pair_coeff 20 20 0.03 2.5 # H6 pair_coeff 21 21 0.03 2.5 # H7 pair_coeff 22 22 0.03 2.5 # H8 pair_coeff 23 23 0.03 2.5 # H9 pair_coeff 24 24 0.17 3.25 # N1 pair_coeff 25 25 0.17 3.25 # N2 pair_coeff 26 26 0.066 3.30 # Z1 pair_coeff 27 27 0.066 3.30 #Z2 pair_coeff 28 28 0.03 2.5 #ZH1 pair_coeff 29 29 0.03 2.5 #ZH2 pair_coeff 30 30 0.03 2.5 #ZH3 pair_coeff 31 31 0.170 3.20 #ZN pair_coeff 32 32 0.0474 3.211 # A1 pair_modify shift yes mix arithmetic bond_coeff 1 386.36 1.39 bond_coeff 2 340 1.08 bond_coeff 3 340 1.08 bond_coeff 4 340 1.09 bond_coeff 5 340 1.09 bond_coeff 6 340 1.09 bond_coeff 7 268 1.529 bond_coeff 8 337 1.466 bond_coeff 9 477 1.315 bond_coeff 10 427 1.378 bond_coeff 11 521 1.341 bond_coeff 13 385.0 1.458 bond_coeff 12 410.0 1.087 bond_coeff 14 650.0 1.157 angle_coeff 1 35 130.9 # C1-C2-H2 angle_coeff 2 70 107.1 # C1-C2-N2 angle_coeff 3 70 108 # C1-N1-C3 angle_coeff 4 70 125.6 # C1-N1-C5 angle_coeff 5 35 130.9 # C2-C1-H1 angle_coeff 6 70 107.1 # C2-C1-N1 angle_coeff 7 70 108 # C2-N2-C3 angle_coeff 8 70 126.8 # C2-N2-C6 angle_coeff 9 70 126.4 # C3-N1-C5 angle_coeff 10 70 126.4 # C3-N2-C6 angle_coeff 11 37.5 109.7 # C4-C6-H3 angle_coeff 12 37.5 109.7 # C4-C6-H4 angle_coeff 13 58.35 112.7 # C4-C6-N2 angle_coeff 14 37.5 110.7 # C6-C4-H10 angle_coeff 15 37.5 110.7 # C6-C4-H8 angle_coeff 16 37.5 110.7 # C6-C4-H9 angle_coeff 17 80 109.47 # F1-B1-F2 angle_coeff 18 80 109.47 # F1-B1-F3 angle_coeff 19 80 109.47 # F1-B1-F4 angle_coeff 20 80 109.47 # F2-B1-F3 angle_coeff 21 80 109.47 # F2-B1-F4 angle_coeff 22 80 109.47 # F3-B1-F4 angle_coeff 23 35 122 # H1-C1-N1 angle_coeff 24 33 107.9 # H10-C4-H8 angle_coeff 25 33 107.9 # H10-C4-H9 angle_coeff 26 35 125.1 # H11-C3-N1 angle_coeff 27 35 125.1 # H11-C3-N2 angle_coeff 28 35 122 # H2-C2-N2 angle_coeff 29 33 108.9 # H3-C6-H4 angle_coeff 30 37.5 107.5 # H3-C6-N2 angle_coeff 31 37.5 107.5 # H4-C6-N2 angle_coeff 32 33 109.8 # H5-C5-H6 angle_coeff 33 33 109.8 # H5-C5-H7 angle_coeff 34 37.5 109.2 # H5-C5-N1 angle_coeff 35 33 109.8 # H6-C5-H7 angle_coeff 36 37.5 109.2 # H6-C5-N1 angle_coeff 37 37.5 109.2 # H7-C5-N1 angle_coeff 38 33 107.9 # H8-C4-H9 angle_coeff 39 70 109.8 # N1-C3-N2 angle_coeff 41 35.03 110.0 angle_coeff 40 80.07 180.0 angle_coeff 42 35.03 110.0 angle_coeff 43 35.03 110.0 angle_coeff 44 30.38 109.5 angle_coeff 44 30.38 109.5 angle_coeff 45 30.38 109.5 angle_coeff 46 30.38 109.5 dihedral_coeff 1 -3.23 0 0 0 # C1-C2-N2-C3 dihedral_coeff 2 -3.23 0 0 0 # C1-C2-N2-C6 dihedral_coeff 3 -3.23 0 0 0 # C2-C1-N1-C3 dihedral_coeff 4 -3.23 0 0 0 # C2-C1-N1-C5 dihedral_coeff 5 -3.23 0 0 0 # C4-C6-N2-C2 dihedral_coeff 6 -3.23 0 0 0 # C4-C6-N2-C3 dihedral_coeff 7 0 0 1.331 0 # H1-C1-C2-H2 dihedral_coeff 8 0 0 0.35 0 # H1-C1-C2-N2 dihedral_coeff 9 0 0 0 0 # H1-C1-N1-C3 dihedral_coeff 10 0 0 0 0 # H1-C1-N1-C5 dihedral_coeff 11 0 0 1.331 0 # H10-C4-C6-H3 dihedral_coeff 12 0 0 1.331 0 # H10-C4-C6-H4 dihedral_coeff 13 0 0 0.35 0 # H10-C4-C6-N2 dihedral_coeff 14 0 0 0 0 # H11-C3-N1-C1 dihedral_coeff 15 0 0 0 0 # H11-C3-N1-C5 dihedral_coeff 16 0 0 0 0 # H11-C3-N2-C2 dihedral_coeff 17 0 0 0 0 # H11-C3-N2-C6 dihedral_coeff 18 0 0 0 0 # H2-C2-N2-C3 dihedral_coeff 19 0 0 0 0 # H2-C2-N2-C6 dihedral_coeff 20 0 0 0 0 # H3-C6-N2-C2 dihedral_coeff 21 0 0 0 0 # H3-C6-N2-C3 dihedral_coeff 22 0 0 0 0 # H4-C6-N2-C2 dihedral_coeff 23 0 0 0 0 # H4-C6-N2-C3 dihedral_coeff 24 0 0 0 0 # H5-C5-N1-C1 dihedral_coeff 25 0 0 0 0 # H5-C5-N1-C3 dihedral_coeff 26 0 0 0 0 # H6-C5-N1-C1 dihedral_coeff 27 0 0 0 0 # H6-C5-N1-C3 dihedral_coeff 28 0 0 0 0 # H7-C5-N1-C1 dihedral_coeff 29 0 0 0 0 # H7-C5-N1-C3 dihedral_coeff 30 0 0 1.331 0 # H8-C4-C6-H3 dihedral_coeff 31 0 0 1.331 0 # H8-C4-C6-H4 dihedral_coeff 32 0 0 0.35 0 # H8-C4-C6-N2 dihedral_coeff 33 0 0 1.331 0 # H9-C4-C6-H3 dihedral_coeff 34 0 0 1.331 0 # H9-C4-C6-H4 dihedral_coeff 35 0 0 0.35 0 # H9-C4-C6-N2 dihedral_coeff 36 0 0 0.35 0 # N1-C1-C2-H2 dihedral_coeff 37 0 0 0 0 # N1-C1-C2-N2 dihedral_coeff 38 0 4.65 0 0 # N1-C3-N2-C2 dihedral_coeff 39 0 4.65 0 0 # N1-C3-N2-C6 dihedral_coeff 40 0 4.65 0 0 # N2-C3-N1-C1 dihedral_coeff 41 0 4.65 0 0 # N2-C3-N1-C5 dihedral_coeff 42 0 1 0 0 dihedral_coeff 43 0 1 0 0 dihedral_coeff 44 0 1 0 0 mass 1 10.811 # B1 mass 2 12.0107 # C1 mass 3 12.0107 # C2 mass 4 12.0107 # C3 mass 5 12.0107 # C4 mass 6 12.0107 # C5 mass 7 12.0107 # C6 mass 8 18.9984 # F1 mass 9 18.9984 # F2 mass 10 18.9984 # F3 mass 11 18.9984 # F4 mass 12 12.0107 # A1 mass 13 1.00794 # H1 mass 14 1.00794 # H10 mass 15 1.00794 # H11 mass 16 1.00794 # H2 mass 17 1.00794 # H3 mass 18 1.00794 # H4 mass 19 1.00794 # H5 mass 20 1.00794 # H6 mass 21 1.00794 # H7 mass 22 1.00794 # H8 mass 23 1.00794 # H9 mass 24 14.0067 # N1 mass 25 14.0067 # N2 mass 26 12.0107 # Z1 mass 27 12.0107 # Z2 mass 28 1.00794 # ZH1 mass 29 1.00794 # ZH2 mass 30 1.00794 # ZH3 mass 31 14.0067 # ZN mass 32 12.0107 # A2 #systems group emim type 2 3 4 5 6 7 13 14 15 16 17 18 19 20 21 22 23 24 25 group bf4 type 1 8 9 10 11 group elec1 type 12 group elec2 type 32 group sol type 1 2 3 4 5 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 group syst type 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 group acn type 1 2 3 4 5 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 group ele type 12 32 group elel type 12 group eler type 32 group elelr molecule 44444 group elell molecule 44445 group elerr molecule 44441 group elerl molecule 44442 #charges set type 1 charge 0.66208 # B1 set type 2 charge -0.192 # C1 set type 3 charge -0.192 # C2 set type 4 charge -0.072 # C3 set type 5 charge -0.28 # C4 set type 6 charge -0.192 # C5 set type 7 charge -0.136 # C6 set type 8 charge -0.36552 # F1 set type 9 charge -0.36552 # F2 set type 10 charge -0.36552 # F3 set type 11 charge -0.36552 # F4 set type 12 charge 0 # A1 set type 13 charge 0.216 # H1 set type 14 charge 0.064 # H10 set type 15 charge 0.168 # H11 set type 16 charge 0.216 # H2 set type 17 charge 0.144 # H3 set type 18 charge 0.144 # H4 set type 19 charge 0.144 # H5 set type 20 charge 0.144 # H6 set type 21 charge 0.144 # H7 set type 22 charge 0.064 # H8 set type 23 charge 0.064 # H9 set type 24 charge 0.176 # N1 set type 25 charge 0.176 # N2 set type 32 charge 0 # A2 set type 26 charge -0.08 # Z1 set type 27 charge 0.46 # Z2 set type 28 charge 0.06 # ZH1 set type 29 charge 0.06 # ZH2 set type 30 charge 0.06 # ZH3 set type 31 charge -0.56 # A2 #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_sol sol temp fix 5 ele setforce 0.0 0.0 0.0 fix 6 sol shake 0.0001 20 0 b 6 7 a 1 2 minimize 1.0e-4 1.0e-6 1000 10000 velocity sol create 298.0 86289 #timestep 1 dump 1 syst custom 1000 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." ############################################################### #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 region region1 block -0.0195 49.9805 -0.02849 49.971501 -51.999001 152.000999 units box compute emim_num_density emim chunk/atom bin/1d z -51.99 0.2 region region1 units box compute bf_num_density bf4 chunk/atom bin/1d z -51.99 0.2 region region1 units box compute acn_num_density acn chunk/atom bin/1d z -51.99 0.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_bf4_density.txt compute overall_density sol chunk/atom bin/1d z -51.999001 0.1 region region1 units box fix dens_overall sol ave/chunk 1000 1 1000 overall_density density/mass ave one file densities.txt dump 3 syst custom 100 dir$r/dump/dump1.lammpstrj id type x y z q dump_modify 3 sort id timestep 0.1 fix 3 sol nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_sol run 10000 unfix 3 write_restart dir$r/finalpart1.rest write_data dir$r/finalpart1.dat #region region1 block 0.036 41.724 -0.0015 42.539 -70.029501 69.970499 units box #compute number_density all chunk/atom bin/1d z center 1 units box #compute charge_density all chunk/atom bin/1d z center 1 units box #fix ave_number all ave/chunk 100 1 100 number_density density/number region myregion timestep 1 fix 3 sol nvt temp 298.0 298.0 100.0 fix_modify 3 temp t_sol run 100000 unfix 3 write_restart dir$r/finalpart2.rest write_data dir$r/finalpart2.dat 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/33kembf.csv fix 4 sol nvt temp 298.0 298.0 100 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_t_sol c_ql c_qr c_qll c_qlr c_qrl c_qrr pe etotal ecoul press spcpu thermo_modify flush yes norm no thermo_modify temp t_sol run 10000000