Dear users and developers:
I would like to simulate the diffusion coefficient of n-pentane at 300K and 1 atm by TraPPE-UA force field with MSD method, but the value I got is higher than the literature (the literature mentioned is 6.14e-9m2/s of pure n-pentane in the same condition, what I got is around 7.4e-9m2/s). Then I trid to repeat the processes given in the literature but stilll got the wrong value. When I saw back in my log file, I found the box shrink when I apply NPT ensemble to control the press in 1 atm, and I got a higher density than the literature (the density in literature is 620.19kg/m3, what I got is about 630kg/m3). Here is my input script and a part of data file:
############input script:
units real
dimension 3
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
read_data fprintf.out
pair_style lj/cut 14.0
pair_coeff 1 1 0.1947 3.75 #epsilon,sigma for (CH3-C5)-(CH3-C5)
pair_coeff 1 2 0.1334 3.85 #epsilon sigma for (CH3-C5)-(CH2-C5)
pair_coeff 2 2 0.09138 3.95 #epsilon sigma for (CH2-C5)-(CH2-C5)
pair_modify tail yes
bond_coeff 1 95.8052 1.54
angle_coeff 1 62.050 114.0
dihedral_coeff 1 0.00 1.4099 -0.27068 3.1425
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
#structure optimization
replicate 2 2 2
thermo 100
thermo_style custom step temp etotal press vol
minimize 1.0e-4 1.0e-6 10000 10000
write_restart restart.%.*
#NVT ensemble: up to 700k
timestep 0.1
fix 1 all nvt temp 0.1 700.0 10.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
run 1000000
unfix 1
keep system in 700k for a while
fix 2 all nvt temp 700.0 700.0 10.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
run 2000000
unfix 2
down the system 500k
fix 3 all nvt temp 700.0 500.0 10.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
run 3000000
unfix 3
#down the system 300k
fix 4 all nvt temp 500.0 300.0 10.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
run 5000000
unfix 4
#keep system in 300k for a while
fix 5 all nvt temp 300.0 300.0 10.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
run 5000000
unfix 5
#NPT
timestep 1.0
thermo 5000
thermo_style custom step temp ke pe etotal press vol
fix presscontrol all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run 5000000
unfix presscontrol
begin control temperature in 300k for 2ns
reset_timestep 0
fix 6 all nvt temp 300.0 300.0 100.0
compute 1 all msd com yes
variable msdC5 equal c_1[4]/6/(step*dt+1.0e-6) #1.0e-6 is for denominator is not zero
thermo 5000
thermo_style custom step temp ke pe etotal press vol c_1[4] v_msdC5
run 2000000
unfix 6
#nve for 5ns
reset_timestep 0
fix 7 all nve
compute 2 all msd com yes
variable msdC52 equal c_2[4]/6/(step*dt+1.0e-6) #1.0e-6 is for denominator is not zero
thermo 5000
thermo_style custom step temp ke pe etotal press vol c_2[4] v_msdC52
run 5000000
unfix 7
EMD for 5ns
reset_timestep 0
fix emd1 all nve
compute 3 all msd com yes
variable msdC53 equal c_3[4]/6/(step*dt+1.0e-6) #1.0e-6 is for denominator is not zero
thermo 5000
thermo_style custom step temp ke pe etotal press vol c_3[4] v_msdC53
run 5000000
#########data file
625 atoms
500 bonds
375 angles
250 dihedrals
2 atom types
1 bond types
1 angle types
1 dihedral types
0.000000 28.921897 xlo xhi
0.000000 28.921897 ylo yhi
0.000000 28.921897 zlo zhi
Masses
1 15.033000
2 14.028000
Atoms
1 1 1 0.000000 0.030030 13.526048 4.639302
… …
625 125 1 0.000000 5.860296 7.088066 9.130796
Bonds
1 1 1 2
… …
500 1 624 625
Angles
1 1 1 2 3
… …
375 1 623 624 625
Dihedrals
1 1 1 2 3 4
… …
250 1 622 623 624 625
I have put the input script (in.C5), data file (fprint.out) and log file in the attachments, please check.
Did the density effect the coefficient a lot? Did I set something wrong in my input script so that I got higher values of density and diffusion coefficient? I have tried to run longer but it does not working well.
The literature is [Antoun S , Saghir M Z , Srinivasan S . An improved molecular dynamics algorithm to study thermodiffusion in binary hydrocarbon mixtures. J. Chem. Phys., 2018, 148 104507.]
Best wishes for you
Xiaoyu
fprintf.out (56.4 KB)
in.C5 (3.29 KB)
log.lammps (649 KB)