#echo screen
units real
variable T equal 298
variable V equal vol
variable dt equal 0.00025
variable p equal 200  #correlation length
variable s equal 10   #sampling interval
variable d equal $p*$s #dump interval 
variable dtPreEquilib equal 0.000005
variable liqx equal 35.000
variable liqy equal 40.000
variable liqz equal 35.000

# unit conversion

variable kB equal 1.3806504e-23 # Boltzmann's constant
variable kCal2j equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal ${kCal2j}*${kCal2j}/${fs2s}/${A2m}

#setup problem

dimension 3
boundary p p f 
atom_style full
read_data gw.data

#create regions for the graphene walls

region 1 block INF INF INF INF INF -6.0 units box           # lower graphene walls 
region 1a block INF INF INF INF INF -2.0 units box        
region 1b block INF INF INF INF -6.0 -2.0 units box

region 2 block INF INF INF INF 36.0 INF  units box           # upper graphene walls
region 2a block INF INF INF INF 36.0 40.0  units box
region 2b block INF INF INF INF 40.0 INF  units box

group lower region 1	# Define group "lower" for lower wall
group lowera region 1a  # Define group "lowera" for the bottom graphene layer of the lower wall
group lowerb region 1b  # Define group "lowerb" for the upper graphene layer of the lower wall

group upper region 2	# Define group "upper" for upper wall
group uppera region 2a	# Define group "uppera" for lower graphene layer of the upper wall
group upperb region 2b	# Define group "upperb" for upper graphene layer of the upper wall


variable liquidvol equal ${liqx}*${liqy}*${liqz} 


# create region of the water in the xz plane for the poiseuille flow 

region 4 block 0.0 35.00 0.0 1.0 0.0 35.0 units box
group flow region 4

#region 5 for the water molecules 

region 5 block INF INF INF INF 0.0 35.0 units box
group liquid region 5


pair_style hybrid lj/charmm/coul/long 10.0 10.5 10.5 &  			  
lj/charmm/coul/charmm 9.0 10.0
pair_coeff 1 1 lj/charmm/coul/long 0.1553 3.166		
pair_coeff 1 2 lj/charmm/coul/long 0.0 2.058
pair_coeff 2 2 lj/charmm/coul/long 0.00774378 0.98
pair_coeff 1 3 lj/charmm/coul/long 0.3135 3.19
pair_coeff 2 3 lj/charmm/coul/long 0.0 3.28
pair_coeff 3 3 lj/charmm/coul/charmm 0.068443 3.407



# bond and angle interactions for the water molecules

bond_style hybrid harmonic
bond_coeff 1 harmonic 1000 1.0

angle_style hybrid harmonic
angle_coeff 1 harmonic 1000 109.47
group spce type 1 2

fix fshake spce shake 0.0001 10 100 b 1 a 1
kspace_style pppm 1e-05
kspace_modify slab 3.0
pair_modify mix arithmetic
neighbor 2.0 bin
neigh_modify delay 1 every 1 check yes
neigh_modify exclude group upper upper
neigh_modify exclude group lower lower


velocity liquid create 300.0 34239


timestep     ${dt} # Sets the timestep for the simulations
thermo	     $d # Prints thermodynamic info every d seconds
velocity     all create $T 102486 mom yes rot yes dist gaussian # Sets the initial temperatures
fix          NVE all nve temp $T $T 10 drag 0.2 # Use the microcanonical ensemble for the simulation


# Sets thermostats for each graphene walls
fix		templwal1a  lowera	temp/rescale	10 298.0 298.0 0.02 1.0
fix		templwal1b  lowerb	temp/rescale	10 298.0 298.0 0.02 1.0
fix		templwal2ba uppera	temp/rescale	10 298.0 298.0 0.02 1.0
fix		templwal2bb upperb	temp/rescale	10 298.0 298.0 0.02 1.0

# Fixes the temperature of the liquid for the equilibriation phase
fix		liquid      liquid      temp/rescale	10 298.0 298.0 0.02 1.0


unfix fshake
minimize 1.0e-4 1.0e-6 100 1000
fix fshake spce shake 0.0001 10 100 b 1 a 1




# Calculates the kinetic energy and prints density, velocity
# and temperature profiles in a file
compute   keLiquid liquid ke/atom
variable  temp atom c_keLiquid/1.5
fix           aveData  liquid ave/spatial 1  200000 200000  z lower   0.005 density/number  vy v_temp units reduced file graph_w.profile



timestep	${dtPreEquilib}
run 10000
dump		1 all atom 50000 dumpfile.lammpstrj

timestep	${dt}
 run          200000 # Equilibriation phase without nanoparticle

reset_timestep 0
# Green-Kubo model for calculating the thermal conductivity of the liquid
compute myKE liquid ke/atom
variable atomTempLiquid atom c_myKE/1.5
compute tempLiquid liquid reduce ave v_atomTempLiquid
compute      myPE liquid pe/atom
compute      myStress liquid stress/atom virial
compute      flux liquid heat/flux myKE myPE myStress
variable     Jx equal c_flux[1]/vol
variable     Jy equal c_flux[2]/vol
variable     Jz equal c_flux[3]/vol
fix          JJ liquid ave/correlate $s $p $d &
             c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running
variable     scale1 equal ${convert}/${kB}/c_tempLiquid/c_tempLiquid/v_liquidvol*$s*${dt}
variable     k11 equal trap(f_JJ[3])*v_scale1
variable     k22 equal trap(f_JJ[4])*v_scale1
variable     k33 equal trap(f_JJ[5])*v_scale1
variable     tempk equal (v_k11+v_k22+v_k33)/3.0 
thermo_style custom step temp v_liquidvol v_tempk
dump		1 all atom 50000 dumpfile.lammpstrj



# variable     scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt}

# calculate viscosity
#reset_timestep 0
#variable     pxy equal pxy
#variable     pxz equal pxz
#variable     pyz equal pyz
#fix          SS liquid ave/correlate $s $p $d &
             v_pxy v_pxz v_pyz type auto file S0St.dat ave running
#variable     scale equal ${convert}/(${kB}*$T)*v_liquidvol*$s*${dt}
#variable     v11 equal trap(f_SS[3])*${scale}
#variable     v22 equal trap(f_SS[4])*${scale}
#variable     v33 equal trap(f_SS[5])*${scale}
#thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
#dump 1 flow custom 2500000 dump2.flow z vy
#dump 2 all atom 50000 flow2.lammpstrj
#fix  ans flow ave/spatial 2 5000 5000000 z center 0.1 z vy density/number file fix2.dat
run 100000 # Run simulation
#variable     v equal (v_v11+v_v22+v_v33)/3.0
#print        "average viscosity: $v [Pa.s] @ $T K"
variable     k equal (v_k11+v_k22+v_k33)/3.0 
print        "average conductivity: $k[W/mK] @ $T K"