# water surface tension 

units		real
atom_style	full
dimension	3
variable        x equal temp
variable        t equal 298

# initialization *********************************************************************************************
read_data       Waterbulk298.data
group           water type 1 2
region          1 block 6.8760546260141009e+00 3.8123945373986558e+01 6.8760546260141009e+00 3.8123945373986558e+01 3.7123945373986558e+01 3.8223945373986558e+01
region          2 block 6.8760546260141009e+00 3.8123945373986558e+01 6.8760546260141009e+00 3.8123945373986558e+01 6.7760546260141009e+00 7.4760546260141009e+00
region          3 block 6.8760546260141009e+00 3.8123945373986558e+01 6.8760546260141009e+00 3.8123945373986558e+01 15 25
delete_atoms	region 1 mol yes
delete_atoms	region 2 mol yes
change_box	all z scale 3
# neighbor **************************************************************************************************

neighbor   	2 multi
neigh_modify	delay 1
#comm_modify    mode multi

# water potential **********************************************************************************
pair_style lj/cut/coul/long 10.0 
pair_coeff 1 1 0 0  #H-H
pair_coeff 1 2 0 0 #H-O
pair_coeff 2 2 0.1553  3.166 #O-O

bond_style harmonic
bond_coeff 1 450.0 1

angle_style harmonic
angle_coeff 1 55 109.47

set             type 1 charge 0.4238
set             type 2 charge -0.8476

kspace_style    ewald 1.0e-6


#fixes***********************************************************************************************

fix		1 water shake 0.0001 30 0 b 1 a 1
fix		2 all nvt temp 298 298 100.0



#compute and variables and thermo******************************************************************************************L
compute		my_pres all pressure thermo_temp
variable	my_gamma equal lz*(c_my_pres[3]-(c_my_pres[1]+c_my_pres[2])/2.0)/2
compute		msd all msd 	   
dump		1 all atom 500 water.xyz
thermo_style	custom step density temp etotal pe ke lz  press v_my_gamma v_dens
thermo		100


#fixes 2***********************************************************************************************
fix             3 all ave/time 50 10  500 c_msd file msd.xls mode vector
fix		4 all ave/time 500 200 200000 v_my_gamma file st.xls mode scalar
fix		5 all ave/time 50 200 20000 v_my_gamma file st2.xls mode scalar


#run ***********************************************************************************************

timestep        1.0
restart		100000	poly.restart
run             600000
write_data      st298.data