units real

atom_style atomic

# substance

variable substance1 string deuterium

variable substance2 string neon

variable molarweight1 equal 4.028209954364

variable molarweight2 equal 20.17969806457545

# thermodynamic state

variable x1 equal 0.2

variable initial_density equal 813.43 # kg / m3 @ p = 8 bar

variable t equal 32

variable natoms equal 4000

variable n1 equal v_x1*v_natoms

variable n2 equal v_natoms-v_n1

variable molarweight_mix equal v_x1*v_molarweight1+(1-v_x1)*v_molarweight2

variable avogadro equal 6.022*10^23

variable initial_density_mol equal v_initial_density/v_molarweight_mix*v_avogadro*10^(-27) # N / V

variable unit_cell_volume equal 1/v_initial_density_mol

variable box_length equal v_unit_cell_volume^(1/3)*v_natoms^(1/3)

variable rc equal 12.0

# setting up the simulation box

region simulation_box block 0.0 ${box_length} 0.0 ${box_length} 0.0 ${box_length}

create_box 2 simulation_box

create_atoms 1 random ${n1} 1234 NULL

create_atoms 2 random ${n2} 4321 NULL

mass 1 ${molarweight1}

mass 2 ${molarweight2}

# pair potentials

pair_style table linear 5000

pair_coeff 1 1 ${substance1}_${t}K.table FH1_${substance1} ${rc}

pair_coeff 2 2 ${substance2}_${t}K.table FH1_${substance2} ${rc}

pair_coeff 1 2 ${substance1}_${substance2}_${t}K.table FH1_${substance1}_${substance2} ${rc}

neighbor 0.3 bin

neigh_modify delay 0 every 1

# energy minimization

print "Running Minimization"

minimize 1.0e-4 1.0e-6 100 1000

# generate temperatures

velocity all create $t 12345

# Define output

thermo_style custom step temp pe ke vol density

print "Running Equilibration"

fix 1 all nvt temp $t $t 500

run 100000

reset_timestep 0

# RDF

compute gr all rdf 500 1 1 1 2 2 1 2 2

fix 1111 all ave/time 100 1 100 c_gr[*] &

file ${substance1}_${substance2}_${t}K.rdf &

mode vector ave running overwrite &

title1 "# g(r) for ${substance1} ${substance2} at rho=${initial_density}, t=${t}" &

title2 "i r g11 c11 g12 c12 g21 c21 g22 c22"

thermo 1000

fix 1 all nvt temp $t $t 500

dump trj all atom 1000 ${substance1}_${substance2}_${t}K.lammpstrj

dump_modify trj sort id

dump xtctrj all xtc 1000 ${substance1}_${substance2}_${t}K.xtc

dump_modify xtctrj sort id

print "Running Production"

run 100000