units real
atom_style full
pair_style lj/cut/coul/long 12.0
kspace_style pppm 1e-5
bond_style harmonic
angle_style harmonic
variable delete_buffer index 1.0 # delete water molecules this distance from the box edge.
variable slab_offset index 5.0
variable temp index 300.0
variable press index 1.0
variable elyte_data index elyte.data
variable slab_left_data index slab.data
variable slab_right_data index slab.data
variable n_additional index 1 # Number of additional atom types to add when reading in the data files.
read_data ${elyte_data} extra/atom/types ${n_additional}
group water type 1 2
# Delete water molecules with bonds over PBCs.
region delete_lo block EDGE EDGE EDGE EDGE $(zlo) $(zlo+v_delete_buffer)
# using random because supports both region and group
delete_atoms random fraction 1.0 yes water delete_lo 41234 mol yes bond yes
region delete_hi block EDGE EDGE EDGE EDGE $(zhi-v_delete_buffer) $(zhi)
# using random because supports both region and group
delete_atoms random fraction 1.0 yes water delete_hi 8080234 mol yes bond yes
# Read the slabs into the system.
read_data ${slab_left_data} shift 0 0 $(zhi+v_slab_offset) add append
read_data ${slab_right_data} shift 0 0 $(zlo-v_slab_offset) add append
write_data pre_min.data
# Minimisation
fix waterbondangles all shake 1.0e-5 20 0 b 1 a 1
fix npt all npt temp $(v_temp) $(v_temp) 100.0 aniso $(v_press) $(v_press) $(1000*dt) couple xy
min_style cg
minimize 1.0e-4 1.0e-6 100 1000
write_data minimized.data
unfix npt
# NVE with slowly increasing interaction strength
variable scale equal ramp(0.001,1.0)
#fix nve all nve
#fix adapt all adapt 10 pair lj/cut/coul/long epsilon * * v_scale scale yes reset yes
#
#fix waterbondangles all shake 1.0e-5 20 0 b 1 a 1
##dump adapt_run all custom 100 adapt.dump id type x y z q
#thermo 100
#run 1000
#unfix adapt
#
#write_data adapt.data
#
#unfix nve
# NVT Run.
fix nvt all nvt temp $(v_temp) $(v_temp) 100.0
fix waterbondangles all shake 1.0e-5 20 0 b 1 a 1
thermo 100
run 1000
unfix nvt
write_data nvt.data
# NPT Run.
fix npt all npt temp $(v_temp) $(v_temp) 100.0 aniso $(v_press) $(v_press) $(1000*dt) couple xy
thermo 100
run 1000
write_data npt.data
group ions type 3 4
variable avogadro equal 6.0221e23
variable molality equal (count(ions)/2/v_avogadro)/(mass(water)*1e-3/v_avogadro)
print "${molality}"
# Quick NPT run to get average density
variable lx equal lx
variable ly equal ly
variable lz equal lz
# Average over the last 5000 steps.
reset_timestep 0 # Because fix only outputs on multiples of n_npt
variable n_npt index 10000
thermo_style custom step temp press lx ly lz pxx pyy pzz density
fix ave_box all ave/time 1 ${n_npt} ${n_npt} v_lx v_ly v_lz
run ${n_npt}
variable scale_x equal f_ave_box[1]/lx
variable scale_y equal f_ave_box[2]/ly
variable scale_z equal f_ave_box[3]/lz
print "${scale_x} ${scale_y} ${scale_z} "
change_box all x scale ${scale_x} y scale ${scale_y} z scale ${scale_z}
write_data averaged_box.data