#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
variable conc equal 3
######################################
units real 
atom_style full
dimension 3
newton on
processors * * 2 
boundary p p f

read_data    input_data/IL1.dat
#read_restart dir$r/input_data/

################################################
############force field########################
############################################

pair_style           lj/cut/coul/long 12
bond_style           harmonic
angle_style          harmonic
kspace_style         pppm/electrode 1.0e-6
kspace_modify        slab 3

fix zwalls all wall/reflect zlo EDGE zhi EDGE units box
# Pair Coeffs
#
# 1  C1
# 2  C2
# 3  C3
# 4  C4
# 5  C5
# 6  C6
# 7  C7
# 8  C8
# 9  F1
# 10  F2
# 11  F3
# 12  F4
# 13  F5
# 14  F6
# 15  H1
# 16  H10
# 17  H11
# 18  H2
# 19  H3
# 20  H4
# 21  H5
# 22  H6
# 23  H7
# 24  H8
# 25  H9
# 26  N1
# 27  N2
# 28  N3
# 29  O1
# 30  O2
# 31  O3
# 32  O4
# 33  S1
# 34  S2
# 35  C9
# 36   C10
# Bond Coeffs
#
# 1  

# Masses

mass 1 12.010700 # C1
mass 2 12.010700 # C2
mass 3 12.010700 # C3
mass 4 12.010700 # C4
mass 5 12.010700 # C5
mass 6 12.010700 # C6
mass 7 12.010700 # C7
mass 8 12.010700 # C8
mass 9 18.998404 # F1
mass 10 18.998404 # F2
mass 11 18.998404 # F3
mass 12 18.998404 # F4
mass 13 18.998404 # F5
mass 14 18.998404 # F6
mass 15 1.007940 # H1
mass 16 1.007940 # H10
mass 17 1.007940 # H11
mass 18 1.007940 # H2
mass 19 1.007940 # H3
mass 20 1.007940 # H4
mass 21 1.007940 # H5
mass 22 1.007940 # H6
mass 23 1.007940 # H7
mass 24 1.007940 # H8
mass 25 1.007940 # H9
mass 26 14.006700 # N1
mass 27 14.006700 # N2
mass 28 14.006700 # N3
mass 29 15.999400 # O1
mass 30 15.999400 # O2
mass 31 15.999400 # O3
mass 32 15.999400 # O4
mass 33 32.064999 # S1
mass 34 32.064999 # S2
mass 35 12.010700 # C9
mass 36 12.010700 # C10

#systems
group c2mim type 1 2 3 4 5 6 26 27 18 19 20 21 22 23 24 25 16 17 15
group ntf2 type 7 8 28 9 10 11 12 13 14 29 30 31 32 33 34 28 
group elec1 type 35 
group elec2 type 36
group sol type 1 2 3 4 5 6 26 27 18 19 20 21 22 23 24 25 7 8 28 9 10 11 12 13 14 16 17 15 29 30 31 32 33 34 28 
group syst type 1 2 3 4 5 6 26 27 18 19 20 21 22 23 24 25 7 8 28 9 10 11 12 13 14 16 17 15 35 29 30 31 32 33 34 28 36
group ele type 35 36

#charges
set type 1 charge -0.192 
set type 2 charge -0.192
set type 3 charge  -0.072
set type 4 charge -0.192
set type 5 charge -0.28
set type 6 charge -0.136
set type 7 charge 0.35
set type 8 charge 0.35
set type 9 charge -0.16
set type 10 charge -0.16
set type 11 charge -0.16
set type 12 charge -0.16 
set type 13 charge -0.16
set type 14 charge -0.16
set type 15 charge 0.216
set type 16 charge 0.064
set type 17 charge 0.168
set type 18 charge 0.064
set type 19 charge 0.144
set type 20 charge 0.144
set type 21 charge 0.144
set type 22 charge 0.144
set type 23 charge 0.144
set type 24 charge 0.064 
set type 25 charge 0.064
set type 26 charge 0.176
set type 27 charge 0.176
set type 28 charge -0.528
set type 29 charge -0.424
set type 30 charge -0.424
set type 31 charge -0.424
set type 32 charge -0.424
set type 33 charge 1.02
set type 34 charge 1.02
set type 35 charge 0.00
set type 36 charge 0.00

#modification
pair_modify 	shift yes mix arithmetic
#bond coeff

bond_coeff      1 450.0 0.9572
bond_coeff      2 450.0 0.9572  
bond_coeff      3 450.0 0.9572
#pair coefficients
pair_coeff	1	1	0.07	3.55
pair_coeff	2	2	0.07	3.55
pair_coeff	3	3	0.07	3.55
pair_coeff	4	4	0.066	3.5
pair_coeff	5	5	0.066	3.5
pair_coeff	6	6	0.066	3.5
pair_coeff	7	7	0.066	3.5
pair_coeff	8	8	0.066	3.5
pair_coeff	9	9	0.053	2.95
pair_coeff	10	10	0.053	2.95
pair_coeff	11	11	0.053	2.95
pair_coeff	12	12	0.053	2.95
pair_coeff	13	13	0.053	2.95
pair_coeff	14	14	0.053	2.95
pair_coeff	15	15	0.03	2.42
pair_coeff	16	16	0.03	2.5
pair_coeff	17	17	0.03	2.42
pair_coeff	18	18	0.03	2.42
pair_coeff	19	19	0.03	2.5
pair_coeff	20	20	0.03	2.5
pair_coeff	21	21	0.03	2.5
pair_coeff	22	22	0.03	2.5
pair_coeff	23	23	0.03	2.5
pair_coeff	24	24	0.17	3.25
pair_coeff	25	25	0.03	2.5
pair_coeff	26	26	0.17	3.25
pair_coeff	27	27	0.17	3.25
pair_coeff	28	28	0.17	3.25
pair_coeff	29	29	0.21	2.96
pair_coeff	30	30	0.21	2.96
pair_coeff	31	31	0.21	2.96
pair_coeff	32	32	0.21	2.96
pair_coeff	33	33	0.25	3.55
pair_coeff	34	34	0.25	3.55
pair_coeff      35      35      0.07    3.55
pair_coeff      36      36      0.07    3.55


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}
minimize             1.0e-4 1.0e-6 10 10000
velocity            sol create 0.1 86289
fix 5 ele setforce 0.0 0.0 0.0
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 100000 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}


timestep 1
fix equi sol nvt temp 294.0 294.0 100.0
run 200000
unfix equi 
fix heating sol nvt temp 294.0 800.0 100.0
unfix heating
run 200000
fix equi sol nvt temp 800.0 800.0 100.0
run 1000000
unfix equi
fix cooling sol nvt temp 800.0 294.0 100.0
run 500000
unfix cooling  
fix equi sol nvt temp 294.0 294.0 100.0
run 200000
unfix equi 

#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 0.1
fix 3 IL nvt temp 294.0 294.0 100.0
run 10000

timestep 0.5
fix 3 IL nvt temp 294.0 294.0 100.0
run 10000

variable            w equal ${a}/2
fix                 e elec1 electrode/conp -${w} 1.979 couple elec2 ${w} symm on etypes on 
fix                 4 sol nvt temp 294.0 294.0 100 
fix_modify          4 temp t_water

timestep 1
fix 3 IL nvt temp 294.0 294.0 100.0
run 10000
#Barostat