# GCMC NVT for CO2 molecular fluid, flexible model

echo		both

# variables available on command line

variable        mu index -8.1
variable        Pressu          equal  300
variable        Precoeff        equal  0.3864
variable	      disp index 0.5
variable        T index 338.0
variable        lbox index 45
variable        spacing index 15.0

variable                R               equal   0.00198722
variable                sysvol          equal   vol
variable                sysmass         equal   mass(all)/6.0221367e+23
variable                sysdensity      equal   v_sysmass/v_sysvol/1.0e-24
variable                coulomb         equal   ecoul+elong
variable                etotal          equal   etotal
variable                pe              equal   pe
variable                ke              equal   ke
variable                evdwl           equal   evdwl
variable                epair           equal   epair
variable                ebond           equal   ebond
variable                eangle          equal   eangle
variable                edihed          equal   edihed
variable                time            equal   step*dt+0.000001



# global model settings

units           real
atom_style      full
boundary        p p p
pair_style      lj/class2/coul/long  9.5 
pair_modify     mix sixthpower tail yes
kspace_style    ewald 0.0001
bond_style      class2
angle_style     class2

# 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 co2mol_gcmcnvt.dat
create_atoms   	0 box mol co2mol 464563 units box
                        
# rigid CO2 TraPPE model 

pair_coeff      1   1  0.060 3.3000
pair_coeff      2   2  0.104 3.5258 
bond_coeff      1     1.1600  1161.3421 -2564.5706  3932.8735  
angle_coeff     1     180.0005    57.1000     0.0000     0.0000 
angle_coeff     1     bb 275.4350     1.1600     1.1600
angle_coeff     1     ba 0.0000     0.0000     1.1600     1.1600

# masses

mass 1 12.011000 
mass 2 15.999400 

# MD settings

group           co2 type 1 2

#-------------------------------------------------------------------------------
# Stage 0: minimization
#
minimize 1.0e-4 1.0e-6 100 1000

compute moltemp all temp
compute_modify moltemp dynamic/dof yes
compute_modify thermo_temp dynamic/dof yes

neighbor        2.0 bin
neigh_modify    every 1 delay 10 check yes
velocity       	all create ${T} 54654
timestep        0.1

fix                     mynvt all nvt temp ${T} ${T} 10
run 500000
unfix mynvt

reset_timestep 0


# NVT thermostat and dynamics

fix                     mynvt all nvt temp ${T} ${T} 10
fix_modify		   mynvt  temp moltemp

# gcmc

variable        tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix             mygcmc co2 gcmc 100 10 0 0 54341 ${T} ${mu} ${disp} pressure ${Pressu} fugacity_coeff ${Precoeff} mol co2mol tfac_insert ${tfac} full_energy


# output

variable	tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable	iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable	dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable	racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
compute_modify  thermo_temp dynamic/dof yes
thermo_style    custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
thermo          1000

fix			instantaneous all ave/time 1000 10 10000 v_time  v_sysdensity c_thermo_temp c_thermo_press v_sysvol v_etotal & 
			v_pe v_ke v_evdwl v_coulomb file Instantaneous_300_gcmc.txt
			
# run

run             50000000