# Created by charmm2lammps v1.9.1 on Mon 13 Dec 2021 12:44:15 AM +0330
# Command: charmm2lammps.pl all27_prot_lipid water1_autopsf

echo            both
units           real
boundary        p p p
neigh_modify    delay 2 every 1
timestep        1
atom_style      full

bond_style      harmonic
angle_style     harmonic

read_data       config.data

set  type 1 charge 0.4238
set  type 2 charge -0.8476

pair_style       lj/charmm/coul/long 8 10 


pair_coeff     1  1     0.0     0.0     #H-H
pair_coeff     2  2    0.1553    3.166     #O-O
pair_coeff     3  3     0.07    3.9848     #C-C

pair_coeff     1  2    0.0     0.0     #H-O
pair_coeff     1  3      0.0     1.9924     #H-C
pair_coeff     2  3      0.104264     3.5754     #O-C


#pair_modify     mix arithmetic
kspace_style    pppm 1e-4

neighbor	      2 bin
neigh_modify     delay 0 every 1 check yes
####################################################


#min_style cg
#minimize 0.0 1.0e-8 1000 10000
#minimize 1.0e-6 1.0e-6 1000 100000

###################################################
group water_molec type 1 2

velocity water_molec create 300.0 12345678 dist uniform  #Apply Nose-Hoover termostate on water molecules

fix NVT water_molec nvt temp 300.0 300.0 100      #Apply Nose-Hoover termostate on water molecules

fix rigid water_molec shake 1.0e-6 500 0 t 1 2 a 1 b 2


group piston id <> 11764  12771 
group cnt id <> 13780  14419
group sheet2 id <> 12772 13779
group sheet34 id <> 14420 16435

group g_all union cnt sheet2 sheet34 

velocity piston   set 0 0 0 units box
velocity g_all     set 0 0 0 units box



fix NOFORCE1 piston setforce 0 0 0
fix NOFORCE2 g_all setforce 0 0 0



dump 1 all custom 10000 dump.lammpstrj id type x y z
run 2000    # 0.1ns 100000 
undump 1
##################

reset_timestep 0

region   cyl_1 cylinder z 25 25 4 59 101 units box
region   cyl_2 cylinder z 25 25 4 159 161 units box

group cyl_1 region cyl_1
group cyl_2 region cyl_2

delete_atoms group cyl_1
delete_atoms group cyl_2


region cntube cylinder z 25 25 10.4752 61 99 units box

############# Applying force to piston ##################

fix push piston nve
fix force piston aveforce 0 0 0.0005 # 1GPa = 0.14393 kcal.mol^-1.A^-3   
#velocity piston   set 0 0 0.00001 units box # 1 m/s = 1e-5 A/fs

##########################################################

region box1 block EDGE EDGE EDGE EDGE 0 60 units box
region box2 block EDGE EDGE EDGE EDGE 100 160 units box
region  cyl_3 cylinder z 25 25 30 61 99 units box

group O type 2


compute 1 O property/atom mass
compute 2 O reduce/region cntube sum c_1
variable oxygcount equal c_2/15.9994
variable OCountcnt equal count(O,cyl_3)

fix count all ave/time 1 10 100 v_oxygcount v_OCountcnt file count_oxygen.profile

thermo 100
thermo_style custom step temp press v_oxygcount v_OCountcnt
################# End Green-Kubo viscosity calculation ###############


dump 2 all custom 1000 dump2.lammpstrj id type x y z
run 1000000 #1ns 1000000