# 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       structure.data

bond_style harmonic	
bond_coeff 2 469.0  1.0     #H-O
bond_coeff 1 450.0  1.418      #C-C

angle_style harmonic
angle_coeff 1  85.0  109.47    #H-O-H

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

kspace_style    pppm 1e-4

neighbor	      2 bin
neigh_modify     exclude type 3 3 
neigh_modify     delay 0 every 1 check yes

#delete_atoms overlap 0.3 all all

group water_molec type 1 2

velocity water_molec create 300.0 12345678 dist uniform 
fix NVT water_molec nvt temp 300.0 300.0 100    

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



#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 100000     # 0.1ns 100000 
#undump 1
##################

reset_timestep 0

region   cyl_1 cylinder z 25 25 6.77 59 101 units box
region   cyl_2 cylinder z 25 25 6.77 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


fix push piston nve
fix force piston aveforce 0 0 0.000364593



region box1 block EDGE EDGE EDGE EDGE 0 60 units box
region box2 block EDGE EDGE EDGE EDGE 100 160 units box
region  cyl_out cylinder z 25 25 6.77 61 99 side in units box
region  cyl_in cylinder z 25 25 3 61 99 side out units box

group O type 2



#################################################################
########## Pressure Tensor calculation #####################

variable p equal 30  #Correlation Length
variable s equal 10    #sample interval
variable d equal $p*$s #dump interval

# convert from LAMMPS real units to SI

variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    atm2Pa equal 101325.0
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    T equal 300
variable    V equal 5756.61
variable    dt equal 1.0
variable    convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}
###########################################################################
#variable pxy equal pxy
#variable pxz equal pxz
#variable pyz equal pyz
#variable pxx equal pxx  
#variable pyy equal pyy  
#variable pzz equal pzz 

group mine1 dynamic O region cyl_out every 5


region cyl_intersect intersect 2 cyl_out cyl_in 
group mine2 dynamic O region cyl_intersect every 5


compute Opress1 mine1 stress/atom NULL

compute Opress2 mine2 stress/atom NULL


variable p4 equal c_Opress1[4]
variable p5 equal c_Opress1[5]
variable p6 equal c_Opress1[6]

fix SS mine1 ave/correlate $s $p $d v_p4 v_p5 v_p6 type auto file S0St.dat ave running

variable     scale equal ${convert}/(${kB}*$T)*$V*$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}
variable     v equal (v_v11+v_v22+v_v33)/3.0


fix visco mine1 ave/time 1 10 30 v_v file visco.profile 

thermo 100
thermo_style custom step temp press v_v
                                    
################# End Pressure Tensor calculation ###############



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