# 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       hybrid lj/charmm/coul/long 9 10 airebo 3.0 


pair_coeff      *  *   airebo CH.airebo NULL NULL C

pair_coeff     1  1   lj/charmm/coul/long   0.0     0.0     #H-H
pair_coeff     2  2   lj/charmm/coul/long   0.1553    3.166     #O-O

pair_coeff     1  2   lj/charmm/coul/long   0.0     0.0     #H-O
pair_coeff     1  3   lj/charmm/coul/long   0.0     0.19924     #H-C
pair_coeff     2  3   lj/charmm/coul/long   0.104264088     1.78224     #O-C

#pair_modify     mix arithmetic
kspace_style    pppm 1e-4


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 m 1.0 a 1

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


velocity piston   set 0 0 0 units box
velocity cnt     set 0 0 0 units box
velocity sheet2   set 0 0 0 units box
velocity sheet34   set 0 0 0 units box


fix NOFORCE1 piston setforce 0 0 0
fix NOFORCE2 cnt setforce 0 0 0
fix NOFORCE3 sheet2 setforce 0 0 0
fix NOFORCE4 sheet34 setforce 0 0 0


thermo 1000
thermo_style custom step temp pe ke  etotal
log logequi.txt 


dump 1 all custom 10000 dump.lammpstrj id type x y z
run 60000    # 0.1ns 100000 
undump 1
##################
reset_timestep 0
unfix NOFORCE1

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 force piston aveforce 0 0 0.0005  # 1GPa = 0.14393 kcal.mol^-1.A^-3     

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

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 110.265779084358 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 1000
thermo_style custom step temp pe ke  etotal
################# End Green-Kubo viscosity calculation ###############


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