# 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 8 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     3  3     lj/charmm/coul/long  0.07    3.4     #C-C

pair_coeff     1  2    lj/charmm/coul/long  0.0     0.0     #H-O
pair_coeff     1  3      lj/charmm/coul/long  0.0     1.7     #H-C
pair_coeff     2  3     lj/charmm/coul/long  0.104264     3.283     #O-C

#pair_modify     mix arithmetic
kspace_style    pppm 1e-4

#delete_atoms overlap 0.3 all all
####################################################


#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 all 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

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

##################################################
###########  Applying periodic electric field ####

region cntube cylinder z 25 25 10.4752 61 99 units box

#variable        E equal 0.05 # V/A
#variable        nu equal 0.000006   #  1GHz = 1.0e-6 cycle/fs
#variable        Ez equal $E*cos(2*PI*${nu}*step)  # Frequency w=2*pi*nu

#fix                 e1 all efield 0 0 v_Ez region cntube 

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

fix push piston nve
fix force piston setforce 0 0 0.0005 # 1GPa = 0.14393 kcal.mol^-1.A^-3   0.000020
#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

#######  NEW ADDED????????????????????????????????????????

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


#################################################################
########## Pressure Tensor #####################


variable p equal 40  #Correlation Length
variable s equal 5    #sample interval
variable d equal $p*$s #dump interval


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 mine dynamic O region cyl_3 every 5

fix P_tensor mine ave/time 1 10 100 v_pxx v_pyy v_pzz v_pxy v_pxz v_pyz file ptensor.profile #ptensor

thermo 100
thermo_style custom step temp press v_pxx v_pyy v_pzz v_pxy v_pxz v_pyz
################# End Green-Kubo viscosity calculation ###############


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