# Ethylene Carbonate 

#---------------Simulation parameters------------
units		real
variable        T equal 323.0 # Temperature in Kelvin
variable        Tinit equal 320.0   # equilibration temperature
variable        dt equal 1.0       # Timestep in femtoseconds
variable        V equal vol      # Volume
variable        s equal 5        # Sampling interval
variable        p equal 400       # Number of correlation steps
variable        d equal ${p}*${s}     # Number of times to repeat the correlation, dump interval

#-------convert from LAMMPS real units to SI----------------------
variable        kB equal 1.3806504e-23  # Boltzmann constant [J/K]
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}

#--------------Setup------------------------------------
atom_style	    full
neighbor        2.0 bin
dimension       3
boundary        p p p
neigh_modify    every 2 delay 10 check yes
timestep        ${dt}
thermo          $d

#-----------------------Define groups---------------------------- 
group           all type 1 2 3 4 5

#----------------------Energy Minimization----------------------
dump              1 all custom 100 overall.lammpstrj id type x y z vx vy vz
min_style         cg
minimize          1e-6 1e-6 1000 1000 
#reset_timestep   0

#--------------------Temperature Initialization-----------------
velocity          all create ${Tinit} 6244325 mom yes rot yes dist gaussian
#fix              1 all shake 0.0001 2000 0 b 1 2 3 4 5  a 1 2 3 4 5 

#--------------------------Relaxation------------------------
compute           EC all temp/com
thermo            1000
thermo_style      custom step c_EC temp density vol press ke pe
thermo_modify     flush yes
fix               2 all  nvt temp $T $T 200.0
#restart          100 restart_file.Relax1  restart_file.Relax2
run               500000 
write_restart     relax.restart
write_data        relax.data
#reset_timestep   0

#-------------------------Equilibrium------------------------
compute           EC1 all temp/com
thermo            1000
thermo_style      custom step c_EC1 temp density vol press ke pe etotal
thermo_modify     flush yes
write_data        equilibrium.data
write_restart     equilibrium.restart

unfix             2
fix               3 all npt temp $T $T 200.0 iso 1 1 2000
run               500000  

unfix              3
velocity          all create $T 6244325 mom yes rot yes dist gaussian
fix                4 all nvt temp $T $T 100.0

#--------------------compute viscosity---------------------------------
# Define variables for stress components
variable        pxy equal pxy
variable        pxz equal pxz
variable        pyz equal pyz

# Compute viscosity using Green-Kubo
fix             SS all ave/correlate $s $p $d & 
                v_pxy v_pxz v_pyz 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

thermo_style    custom step temp  press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run            1000000
variable     ndens equal count(all)/vol
print        "average viscosity: $v [Pa.s] @ $mytemp K, ${ndens} atoms/A^3"
      
#--------------------Output Simulation Settings-------------------------
dump               mydump all custom 100 EC.lammpstrj id type x y z 
dump               Dump  all xyz 100 dump.xyz 
dump_modify        Dump  element  c2 oz cz oo hc 
thermo             1000
thermo_style       custom  step temp ke pe etotal density press vol v_v
thermo_modify        flush yes 

write_data my_system.data pair ij