###################################################################
# Input file for simulating (rigid) CO2 using the EPM-2 force field 
#

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

units           real
atom_style      full
boundary        p p p
pair_style      lj/cut/coul/long  14
kspace_style    ewald 1e-6                  # ewald is slower than ppm
pair_modify     mix arithmetic tail yes
bond_style      harmonic
angle_style     harmonic

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

variable        Temp		equal		323
variable        Pres		equal		98.69
variable        spacing     index       5.0         # [Ang] - lattice spaing    

variable        Ninve       equal       500         # Nve initialize
variable        Npnve       equal       1000        # Nve production

variable        Ninvt       equal       500         # Nvt initialize
variable        Npnvt       equal       1000        # Nvt production

variable        Ninpt       equal       500         # Npt initialize
variable        Npnpt       equal       1000        # Npt production

variable        Nf          equal       50          # rate of outputting to the log file 

variable        boxEdgex    index       15.0        # [ang] - box dimensions along x
variable        boxEdgey    index       15.0        # [ang] - box dimensions along y
variable        boxEdgez    index       15.0        # [ang] - box dimensions along z

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

lattice 	sc ${spacing}
region          box block 0 ${boxEdgex} 0 ${boxEdgey} 0 ${boxEdgez} 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

#-------------
# TraPPE model
#-------------

pair_coeff      1   1   0.0536  2.800       # O=[C]=O
pair_coeff      2   2   0.1569  3.050       # [O]=C=O

bond_coeff      1       0       1.16        # [O=C]=O
angle_coeff     1       0       180         # [O=C=O]

#-------------------
# Masses and charges 
#-------------------

mass 1 12.0107      # O=[C]=O
mass 2 15.9994      # [O]=C=O

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

group           co2 type 1 2
neighbor        3.0 bin
neigh_modify    every 1 delay 0 check yes
run_style       verlet

#----------------
# NVE simulations 
#----------------

fix             myrigidnve co2 rigid/nve molecule 

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

dump trjec all atom ${Nf} movie.lammpstrj 

compute 1 co2   bond/local  dist
compute 2 co2   angle/local theta

compute 3 all   reduce ave c_1
compute 4 all   reduce ave c_2

thermo_style custom step pe ke etotal press c_3 c_4 
thermo ${Nf}

variable Temperature    equal   temp
variable Density        equal   density
variable Pressure       equal   press
variable volume         equal   vol

fix     1  all ave/time 4 25 100 v_Temperature v_Density v_Pressure v_volume file thermo.out format %10.6g

#--------------------
# 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.0001
run         ${Ninve} 

timestep    0.001
run         ${Ninve} 

timestep    0.01
run         ${Ninve} 

timestep    0.1
run         ${Ninve} 

timestep    0.2
run         ${Ninve} 

timestep    0.25
run         ${Ninve} 

timestep    0.5
run         ${Ninve} 

timestep    1.0
run         ${Ninve} 

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

timestep    1.0
run         ${Npnve}

write_restart restart.equil.nve

#----
# NVT 
#----

unfix   myrigidnve
fix     myrigidnvt co2 rigid/nvt molecule temp ${Temp} ${Temp} 100         

run ${Ninvt}

run ${Npnvt}

write_restart restart.equil.nvt

#----
# NPT 
#----

unfix   myrigidnvt
fix     myrigidnpt co2 rigid/npt molecule temp ${Temp} ${Temp} 100 iso ${Pres} ${Pres} 1000	

run ${Ninpt}

run ${Npnpt}

write_restart restart.equil.npt