variable 	a loop 1 
variable 	l index 300
log 		log2.${l}
echo         both



#***************VARIABLES**************************#
variable	dt equal 0.0005  #given 0.5 fs 
variable   	T equal 300 

variable  Nevery  equal 10
variable  Nrepeat equal  10
variable  Nfreq equal  100
variable Nthermo equal 100
variable Nstrelax equal 1000
variable Nndrelax equal 5000
variable pressdiff equal 2e8
variable confactor equal 0.6242e9

variable ly equal 40.3
variable lz equal 46.5

dimension       3
units           metal
boundary        p p p
atom_style      full
pair_style 		hybrid lj/charmm/coul/long 9.0 10.0 12.0 lj/cut 10.0

kspace_style 	pppm 1.0e-4
 
#special_bonds 	charmm                                               #added
bond_style   	harmonic
angle_style  	harmonic
pair_modify 	tail yes #table 0

neighbor    	2.0  bin                                              #added
neigh_modify  	every 1 delay 0 check yes #one 20000 page 1000000
 
#******************SET UP*********************
region          mybox block 0 3.1 0 3.1 0 3.1 units box		
create_box      3 mybox bond/types 2 angle/types 2 dihedral/types 1 &
improper/types 1 extra/bond/per/atom 5 extra/angle/per/atom 4 &
extra/dihedral/per/atom 12  extra/improper/per/atom 6 
read_data     "data.singleTIP4P" add append
replicate        43 13 15


region revone block  0     40  0 40.3 0 46.5 units box
region mid block     40    85.36 0 40.3 0 46.5 units box
region revtwo block  85.36 125.36 0 40.3 0 46.5 units box 

group watregone region revone
group watregtwo region revtwo
group watreg union watregone watregtwo
group nonwatreg subtract all watreg
 
delete_atoms group nonwatreg mol yes bond yes
 
 
read_data "graphin" add append shift 42.5 0.7 18.2 group graphinLow        #actual move is 42.5-1.84=40.66                  
read_data "graphin" add append shift 42.5 0.7 28.2 group graphinUp          #create a gap of 0.66 ang
 

#***********PLACING "graphenebound" in LEFT-UP************************

read_data "graphenebound" add append shift 1.84 0.7 0.0  group graboundLftUp
displace_atoms graboundLftUp rotate 0.0 0.0 0.0 1.0 0.0 0.0 90.0 units box
displace_atoms graboundLftUp rotate 0.0 0.0 0.0 0.0 0.0 1.0 90.0 units box
displace_atoms graboundLftUp move 40.7 0.0 28.2 units box 
 
#*********PLACING "graphenebound" in LEFT-DOWN***************************

read_data "graphenebound" add append shift 1.84 0.7 0.0  group graboundLftDown
displace_atoms graboundLftDown rotate 0.0 0.0 0.0 1.0 0.0 0.0 90.0 units box
displace_atoms graboundLftDown rotate 0.0 0.0 0.0 0.0 0.0 1.0 90.0 units box
displace_atoms graboundLftDown move 40.7 0.0 -3.07 units box 
 
#*********PLACING "graphenebound" in RIGHT-UP***************************

read_data "graphenebound" add append shift 1.84 0.7 0.0  group graboundRghtUp
displace_atoms graboundRghtUp rotate 0.0 0.0 0.0 1.0 0.0 0.0 90.0 units box
displace_atoms graboundRghtUp rotate 0.0 0.0 0.0 0.0 0.0 1.0 90.0 units box
displace_atoms graboundRghtUp move 84.8 0.0 28.2 units box 

#*********PLACING "graphenebound" in RIGHT-DOWN***************************

read_data "graphenebound" add append shift 1.84 0.7 0.0  group graboundRghtDown
displace_atoms graboundRghtDown rotate 0.0 0.0 0.0 1.0 0.0 0.0 90.0 units box
displace_atoms graboundRghtDown rotate 0.0 0.0 0.0 0.0 0.0 1.0 90.0 units box
displace_atoms graboundRghtDown move 84.8 0.0 -3.07 units box

group watregone delete
group watregtwo delete
group watreg delete
group nonwatreg delete


#****************MASS*********************************************
mass            3 1.008                                    #H
mass            2 15.9999990                               #O
mass            1 12                                       #C

#***************INTERACTION*************************

pair_coeff  1 1 lj/cut 0.006956 3.4745050026 8.5                              #C-C
pair_coeff  2 2 lj/charmm/coul/long 0.00442314 3.188                          #O-O
pair_coeff  3 3 lj/charmm/coul/long 0.0 0.0                            		  #H-H
pair_coeff  2 3 lj/charmm/coul/long 0.0 0.0                                   #H-O
pair_coeff 	1 2 lj/charmm/coul/long 0.004063 3.19                             #C-O 
pair_coeff  1 3 lj/charmm/coul/long 0.0 3.28                                  #C-H  

bond_coeff   1  43.36  1.5260                                                  # graphene  bond type 1
angle_coeff  1  43.36  110.50                                                  # graphene  angle type 1
bond_coeff   2  19.51  0.9572                                                  # water     bond type 2
angle_coeff  2  2.39   104.52                                                   # water     angle type 2

group  		watadd type 2
group 		spce type 2 3
 
group  		midynagp dynamic watadd region mid every 10
group  		reservL dynamic watadd region revone every 10
group  		reservR dynamic watadd region revtwo every 10
group  		graphene union graphinLow graphinUp
group  		Leftboundgraph union graboundLftUp  graboundLftDown  
group  		Rightboundgraph union  graboundRghtUp graboundRghtDown

timestep ${dt}

neigh_modify  exclude group  graphene graphene                                   
neigh_modify  exclude group  Leftboundgraph Leftboundgraph
neigh_modify  exclude group  Rightboundgraph Rightboundgraph
neigh_modify  exclude group  graphene Leftboundgraph
neigh_modify  exclude group  graphene Rightboundgraph

fix rigid1 graphene setforce 0.0 0.0 0.0
fix rigid2 Leftboundgraph setforce 0.0 0.0 0.0 
fix rigid3 Rightboundgraph setforce 0.0 0.0 0.0 

 
#*****INITIAL VELOCITY**************
velocity all create ${T} 102486 dist gaussian loop geom
 
 #************MINIMIZE*************
thermo 200
thermo_style custom step temp vol press density
thermo_modify lost ignore flush yes
dump 1 all atom 200 dump_min.lammpstrj
minimize 1.0E-9 1.0E-9 10000 10000                                             
reset_timestep 	0

fix   constrain spce shake 1.0e-4 100 0 b 2 a 2
#fix  removeMomentum spce momentum 1 linear 1 1 1 angular
fix   1 spce nvt temp 300 300 0.05                          #100 timestweps  tdamp remember
  
#********TEMPERATURE************************

compute temp_L reservL temp/com
compute temp_R reservR temp/com
compute_modify temp_L dynamic yes
compute_modify temp_R dynamic yes

fix TempAveL reservL ave/time ${Nevery} ${Nrepeat} ${Nfreq} c_temp_L          
fix TempAveR reservR ave/time ${Nevery} ${Nrepeat} ${Nfreq} c_temp_R 

#*************FIRST_RELAXATION************* 
unfix 1
fix 2 spce nve 
thermo ${Nfreq}
thermo_style 	custom step pe ke etotal press c_temp_L c_temp_R f_TempAveL f_TempAveR
thermo_modify   temp temp_L temp temp_R 
dump strelax all atom ${Nthermo} first_relax.lammpstrj
run ${Nstrelax}
undump strelax

#dump 			aftfix2 all atom 200 addforcen.lammpstrj
#run 10000


#*************SECOND_RELAXATION************* 
#********APPLYING ADDFORCE******************
variable 		natoms equal count(reservL,revone)
variable 		n equal ${natoms}/3                                         #varible n is number of oxygen-atoms
variable 		f equal ${pressdiff}*${ly}*${lz}*${confactor}/${n}                                                                                                        #d4.46*10^(-12) n/m^2 only this amount of force                                                                                                         
compute 		peatom reservL pe/atom 
variable 		peatomL atom c_peatom										#addforce in "Nevery" timesteps
fix 			2 watadd addforce v_f 0.0 0.0 energy v_peatomL region revone
dump ndrelax all atom ${Nthermo} second_relax.lammpstrj
run ${Nndrelax}