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

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


angle_coeff 1  85.0  109.47    #H-O-H


set  type 1 charge 0.4238
set  type 2 charge -0.8476





pair_style     hybrid lj/charmm/coul/long 8 10 tersoff
pair_coeff       *  * tersoff SiC_real.tersoff NULL NULL C        #C-C

pair_coeff       1  1      lj/charmm/coul/long      0.046   0.400013524445              # H-H
pair_coeff       2  2      lj/charmm/coul/long      0.1521    3.15057422683             # O-O
pair_coeff       1  2      lj/charmm/coul/long      0.0836    1.7753                    # H-O
pair_coeff       1  3      lj/charmm/coul/long      0.0567    1.9750                    # H-C
pair_coeff       2  3      lj/charmm/coul/long      0.1032    3.503                     # O-C

kspace_style    pppm 1e-6


neighbor	      2 bin
neigh_modify     delay 0 every 1 check yes


min_style cg
minimize 1.0e-4 1.0e-6 100 1000


group solvent type 1 2
group cnt type 3


velocity solvent create 298.0 12345678 dist uniform  #Apply Nose-Hoover termostate on water molecules



delete_atoms overlap 0.3 solvent cnt mol yes


compute mytemp solvent temp
fix 3 solvent nvt temp 298.0 298.0 100
fix_modify 3 temp mytemp


fix push cnt nve
velocity cnt     set 0 0 0 units box
fix NOFORCE1 cnt setforce 0 0 0

region   trash cylinder z 20 20 7.54 1 41 units box side out #7.84
delete_atoms region trash bond yes

thermo 100 
thermo_style custom step temp c_mytemp press pe ke etotal
log logequi.txt  
thermo_modify temp mytemp

dump 1 all custom 1000 dump.lammpstrj id type x y z 
write_data output.data
run 100000 # 100000
undump 1
reset_timestep 0

############################################## Pressure Tensor#######################
group solvent type 1 2
region   cntube cylinder z 20 20 6.77 1 41 units box


compute vis solvent stress/atom mytemp
compute visco solvent reduce sum c_vis[1] c_vis[2] c_vis[3] &
                                  c_vis[4] c_vis[5] c_vis[6]

variable pxxcyl equal -(c_visco[1])/(3*5759.532077)
variable pyycyl equal -(c_visco[2])/(3*5759.532077)
variable pzzcyl equal -(c_visco[3])/(3*5759.532077)
variable pxycyl equal -(c_visco[4])/(3*5759.532077)
variable pxzcyl equal -(c_visco[5])/(3*5759.532077)
variable pyzcyl equal -(c_visco[6])/(3*5759.532077)

group oxyg type 2
variable Countcyl equal count(oxyg,cntube)
thermo 100
thermo_style custom step press temp v_pxxcyl v_pyycyl v_pzzcyl v_pxycyl &
                                    v_pxzcyl v_pyzcyl v_Countcyl
log logequi2.txt 
#thermo_modify press vis
thermo_modify temp mytemp

dump 2 all custom 10000 dump2.lammpstrj id type x y z 
run 50000 # 50000