echo			both

#package omp 4
#suffix omp

dimension        3
boundary	   p p p
units          real
atom_style     full
timestep        1.0             # Timestep/fs
neighbor		        2.0        bin
neigh_modify		    delay 0     every 10     check yes 

##########################################
############ GCMC setting #############
##########################################


variable      mu        equal     		-7.2615024






   
##########################################
############ General setting #############
##########################################

variable        Pdamp           equal     100
variable        seed            equal     649345
variable        rcut            equal      14           # Radius cutoff/A
variable        kspace          equal     0.0001
variable        etol            equal     1e-25
variable        ftol            equal     1e-25
variable        maxiter         equal     100000
variable        maxeval         equal     100000
variable        xLo             equal     xlo
variable        xHi             equal     xhi
variable        yLo             equal     ylo
variable        yHi             equal     yhi
variable        zLo             equal     zlo        
variable        zHi             equal     zhi
variable        dump            equal     100000
variable        thermo          equal     500



##########################################
############ GCMC setting #############
##########################################

#variable        mu              equal     -7.37
#variable       press           equal     100
#variable       fugacity        equal     0.8263 
variable        GCMCNevery      equal     1000
variable        GCMCexchange    equal     100
variable        GCMCmove        equal     0
variable	    disp            equal     0.8
variable        temp            equal     333.15          # Tempereture/K
variable        Tdamp           equal     100
##########################################
############             #############
##########################################




pair_style          lj/cut     ${rcut}                #cutoff lj = cutoff coulomb       
pair_modify         mix                 arithmetic             #Lorentz–Berthelot combining rules     
pair_modify         tail                no            
#kspace_style        pppm                ${kspace}      
#kspace_modify       slab                3.0               
bond_style          harmonic
angle_style         harmonic
dihedral_style      harmonic
improper_style      cvff

read_data              kerogen.data    &  
extra/atom/types  1                &
extra/bond/types   1               &
extra/bond/per/atom 2             &
extra/special/per/atom 2


molecule        C2H6mol C2H6.txt   toff  16  boff 32


pair_coeff      17   17   0.194726   3.75      #ch3  
bond_coeff      33       0       1.54         #ch3-ch3

mass            17        15.035              #ch3


group          kerogen                type     1:16                 
group          C2H6                    type     17



fix          rigidstructure      kerogen    setforce     0    0     0

velocity     kerogen   set    0   0    0


neigh_modify     exclude      group       kerogen   kerogen

delete_bonds kerogen multi any




variable 	typeC2H6 atom "type==17"
group 	        C2H6molecule dynamic C2H6 var typeC2H6
variable        numC2H6      equal     count(C2H6molecule)/2



#-------------------------------------------------------------------------------
# Stage 10: GCMC/NVT ensemble/10 ns/1 fs/  303.15 K- P=100 atm
#-------------------------------------------------------------------------------

fix           myrigidnvt     C2H6    rigid/nvt/small molecule    temp ${temp} ${temp} ${Tdamp}  mol C2H6mol 
fix_modify	myrigidnvt dynamic/dof yes


fix     mygcmc    C2H6    gcmc  ${GCMCNevery}    ${GCMCexchange}   ${GCMCmove}   0   ${seed}  &
${temp}    ${mu}    ${disp}  mol   C2H6mol   group    C2H6   rigid myrigidnvt   tfac_insert  1      #  full_energy      


compute_modify thermo_temp dynamic/dof yes


variable	iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable	dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)

thermo_style     custom    step     v_numC2H6         v_iacc  v_dacc       pe   ke         
thermo ${thermo}


run   5000000              #5 ns