# Created by charmm2lammps v1.9.1 on Thu 03 Feb 2022 06:57:18 PM +0330
# Command: charmm2lammps.pl all27_prot_lipid solvate

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

#pair_style lj/cut/tip4p/long 2 1 1 1 0.1546 13.0
pair_style      lj/charmm/coul/long 8 10

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0	
pair_coeff 2 2 0.1553    3.166     #O-O

kspace_style    pppm 1e-4

fix NVT all nvt temp 300.0 300.0 100

fix rigid all shake 1.0e-6 500 0 m 1.0 a 1

velocity all create 300.0 12345678 dist uniform


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

dump 1 all custom 1000 dump1.lammpstrj id type x y z
run 2000    # 8000  1ns 1000000 
reset_timestep 0

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

variable p equal 40  #Correlation Length
variable s equal 5    #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 convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}


variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
variable pxx equal pxx  #pxx
variable pyy equal pyy  #pyy
variable pzz equal pzz  #pzz


fix SS all ave/correlate $s $p $d &
v_pxy v_pxz v_pyz type auto file S0St.dat ave running

variable T equal 300 
variable V equal vol
variable dt equal 1
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}

thermo 200
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
log logequi2.txt 


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

group O type 2
region box block EDGE EDGE EDGE EDGE EDGE EDGE units box
group mine1 dynamic O region box every 5
compute Opress1 mine1 stress/atom NULL
compute ptensor1 mine1 reduce sum c_Opress1[1] c_Opress1[2] c_Opress1[3] &
                                  c_Opress1[4] c_Opress1[5] c_Opress1[6]

fix P_tensor1 mine1 ave/time 1 5 40 c_ptensor1[1] c_ptensor1[2] c_ptensor1[3] &
                                     c_ptensor1[4] c_ptensor1[5] c_ptensor1[6] file ptensor1.profile #ptensor1

variable v equal (v_v11+v_v22+v_v33)/3.0
variable ndens equal count(all)/vol
print "average viscosity : $v [Pa.s] @ $T K, ${ndens}/A^3}" # ,${vol} A^3 

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