########################################################################################
# Input file for simulating (rigid) CO2 molecules
# Reference : 
# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
# mixtures containing alkanes, carbon dioxide and
# nitrogen AIChE J., 47,1676-1682 (2001)].

#------------------------
# Global model settings 
#------------------------

units           real
atom_style      full
boundary        p p p
pair_style      lj/cut/coul/long  14 
pair_modify     mix arithmetic tail yes
kspace_style    ewald 0.0001
bond_style      harmonic
angle_style     harmonic

#------------------------
# Constants and variables
#------------------------

variable        Temp		equal		338.0       # in K
variable        Pres		equal		1		    # in atm
variable        spacing     index       15.0         # [Ang] - lattice spaing    

variable        Ninvt       equal       10000       # Initialize the NVT ensemble (very small timesteps)
variable        Npnvt       equal       500000      # Production the NVT ensemble
variable        Nf          equal       500         # rate of outputting to the log file 

variable        lbox        index       30.0        # [ang] - box dimensions in angstroms 

#---------------------------------------------
# Box, start molecules on simple cubic lattice 
#---------------------------------------------

lattice 	sc ${spacing}
region          box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
create_box      2 box                       &
                bond/types 1                &
                angle/types 1               &
                extra/bond/per/atom 2       &
                extra/angle/per/atom 1      &
                extra/special/per/atom 2
molecule        co2mol CO2.txt
create_atoms   	0 box mol co2mol 464563 units box

#-----------------------
# rigid CO2 TraPPE model
#-----------------------

pair_coeff      1   1  0.053649   2.8
pair_coeff      2   2  0.156973   3.05
bond_coeff      1       0       1.16
angle_coeff     1       0       180

#-------
# Masses
#-------

mass 1 12.0107 
mass 2 15.9994 

#------------
# MD settings 
#------------

group           co2 type 1 2
neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes
timestep        1.0
run_style       verlet

#------------------------
# NVT simulations - rigid
#------------------------

fix             myrigidnvt co2 rigid/nvt/small molecule temp ${Temp} ${Temp} 10
# fix             myrigidnvt co2 rigid/nvt molecule temp ${Temp} ${Temp} 10

#-------
# Output
#-------

dump trjec all custom 1 movie1.lammpstrj id type x y z vx vy vz 

compute 1 co2 ke            
compute 2 co2 ke/rigid      myrigidnvt      # (scalar) sum of the translational kinetic energy of rigid bodies
compute 3 co2 erotate/rigid myrigidnvt      # (scalar) sum of the rotational    kinetic energy of rigid bodies
compute 4 co2 temp

thermo_style custom step temp pe c_1 c_2 c_3 c_4 ke press vol density
thermo ${Nf}

thermo_modify temp 4

#--------------------
# Run (Initialization)
#--------------------

# initialize to imposed temperature
velocity    all create ${Temp} 54654 
run         0
velocity    all scale ${Temp}
velocity    all zero linear         # set box momentum to zero

timestep    0.01
run         ${Ninvt} 

timestep    0.1
run         ${Ninvt} 

timestep    0.15
run         ${Ninvt} 

timestep    0.2
run         ${Ninvt} 

#-----------------
# Run (Production)
#-----------------

fix  meanTemp  all  ave/time 4 250 1000 c_4 file meanTemp.dat       # calculate a mean temperature

timestep 0.25
run ${Npnvt}