variable disp index 0.5 variable temp index 333.0 variable p1MPa equal 2 variable p2MPa equal 1 variable p1 equal 2*9.869 variable p2 equal 1*9.869 variable lbox index 30.0 variable lbox1 index 10.0 variable spacing index 10.0 variable NAvogad equal 6.02e+23 variable MPatoPa equal 1.0e+6 variable A3tom3 equal 1.0e-30 variable Fvol equal 3000 variable unit2STD equal ${MPatoPa}*${A3tom3} variable Rgasconst equal 8.3144598 variable nmolegas1 equal (${unit2STD}*${p1MPa}*${Fvol})/(${Rgasconst}*${temp}) variable nidealgas1 equal floor(${nmolegas1}*${NAvogad}) variable nmolegas2 equal (${unit2STD}*${p2MPa}*${Fvol})/(${Rgasconst}*${temp}) variable nidealgas2 equal floor(${nmolegas2}*${NAvogad}) # global model settings units real atom_style full boundary f p p pair_style lj/cut 14 #pair_modify mix arithmetic tail yes # box, start molecules on simple cubic lattice lattice sc ${spacing} region box block 0 ${lbox} 0 ${lbox1} 0 ${lbox1} region 1 block 0 ${lbox1} 0 ${lbox1} 0 ${lbox1} region 2 block 20 ${lbox} 0 ${lbox1} 0 ${lbox1} region 3 block 10 20 0 ${lbox1} 0 ${lbox1} create_box 2 box create_atoms 1 random 1 464563 1 create_atoms 2 random 1 999999 2 pair_coeff 1 1 0.2935 3.73 pair_coeff 1 2 0.2935 3.73 pair_coeff 2 2 0.2935 3.73 mass 1 16.042 mass 2 16.042 group 1 region 1 group 2 region 2 neighbor 2.0 bin neigh_modify every 1 delay 0 check yes velocity all create ${temp} 54654 timestep 1 minimize 1.0e-4 1.0e-6 100 1000 reset_timestep 0 fix gcmc1 1 gcmc 1000 100 0 1 99999 ${temp} -0.08 ${disp} region 1 pressure ${p1} fugacity_coeff 0.93128 group 1 grouptype 1 1 tfac_insert 1 #overlap_cutoff 2.5 fix gcmc2 2 gcmc 1000 100 0 2 99999 ${temp} -0.08 ${disp} region 2 pressure ${p2} fugacity_coeff 0.9439 group 2 grouptype 2 2 tfac_insert 1 #overlap_cutoff 2.5 compute t1 1 temp/region 1 compute_modify t1 dynamic/dof yes compute t2 2 temp/region 2 compute_modify t2 dynamic/dof yes variable 1 atom "type==1" variable 2 atom "type==2" group 3 dynamic all region 1 every 1 group 4 dynamic all region 2 every 1 compute p1 3 stress/atom t1 compute p1rd3 3 reduce sum c_p1[1] c_p1[2] c_p1[3] variable press1 equal (-(c_p1rd3[1]+c_p1rd3[2]+c_p1rd3[3])/(3*(vol/3)))*0.101325 compute p2 4 stress/atom t2 compute p2rd3 4 reduce sum c_p2[1] c_p2[2] c_p2[3] variable press2 equal (-(c_p2rd3[1]+c_p2rd3[2]+c_p2rd3[3])/(3*(vol/3)))*0.101325 fix p1 3 ave/time 1 100 100 v_press1 ave running start 2000000 title1 "Average Press[MPa]" file t48p1.txt fix p2 4 ave/time 1 100 100 v_press2 ave running start 2000000 title1 "Average Press[MPa]" file t48p2.txt variable p1ave equal f_p1 variable p2ave equal f_p2 variable 3 equal count(3) variable 4 equal count(4) thermo_style custom step temp press pe ke density atoms v_3 v_4 c_t1 c_t2 v_press1 v_press2 v_p1ave v_p2ave thermo 1000 dump 1 all atom 100000 21.dump run 300000 unfix gcmc1 unfix gcmc2 undump 1 unfix p1 unfix p2 reset_timestep 0 fix gcmc1 1 gcmc 1 10 0 1 99999 ${temp} -0.08 ${disp} region 1 pressure ${p1} fugacity_coeff 0.93128 group 1 grouptype 1 1 tfac_insert 1 fix gcmc2 2 gcmc 1 10 0 2 99999 ${temp} -0.08 ${disp} region 2 pressure ${p2} fugacity_coeff 0.9439 group 2 grouptype 2 2 tfac_insert 1 fix nvt1 all nvt temp 333 333 100 compute t3 1 temp/region 1 compute_modify t3 dynamic/dof yes compute t4 2 temp/region 2 compute_modify t4 dynamic/dof yes compute p3 3 stress/atom t3 compute p3rd3 3 reduce sum c_p3[1] c_p3[2] c_p3[3] variable press3 equal (-(c_p3rd3[1]+c_p3rd3[2]+c_p3rd3[3])/(3*(vol/3)))*0.101325 compute p4 4 stress/atom t4 compute p4rd3 4 reduce sum c_p4[1] c_p4[2] c_p4[3] variable press4 equal (-(c_p4rd3[1]+c_p4rd3[2]+c_p4rd3[3])/(3*(vol/3)))*0.101325 fix p5 3 ave/time 1 100 100 v_press3 ave running start 3000000 title1 "Average Press[MPa]" file t48p5.txt fix p6 4 ave/time 1 100 100 v_press4 ave running start 3000000 title1 "Average Press[MPa]" file t48p6.txt variable p5ave equal f_p5 variable p6ave equal f_p6 thermo_style custom step temp press pe ke density atoms v_3 v_4 v_press3 v_press4 thermo 100 thermo_modify lost ignore dump 3 all atom 1000 32.dump run 1000000