Hello everyone!
I am trying to use RNEMD to determine phonon based thermal conductivity for Copper. After some literature review, I found that the thermal conductivity for Copper at 500K and 300K are ~5 Wm/K and ~10 Wm/K respectively. The dimension used in those studies are roughly similar to what I have used:
36 A x 36 A x 84 A
However, after having done multiple simulations using the Muller Plathe model and following the example script on LAMMPS examples list, the thermal conductivity that I manage for 300K is always round about ~0.6 - 0.7 W/mK. In the example script LJ units were used and I used metal units, so during the final calculations, I convert units of kappa = dQ x dZ/dTemp accordingly.
I initially equilibrated the system containing 10776 atoms at 300K for 5 ns. I used fix npt and the following command:
fix 1 all npt temp 300 300 0.1 iso 1.01325 1.01325 1 #1.01325 bar = 1 atm #Pdamp = 1000*timesteps, but entered as time units, so here 1ps
timestep 0.001
Could it have been that the use of a copper system thermalized at 300K with 1 atm pressure is having an effect on the thermal conductivity results when I am using this thermalized data file for muller plathe model?
I am using the following code for kinetic energy swap of atoms and induction of temperature gradient to eventually calculate thermal conductivity: (total 15.5ns run)
units metal
atom_style atomic
boundary p p p
read_data Data10DT300.R1
pair_style eam #Potential
pair_coeff 1 1 ./Cu_u3.eam
velocity all scale 300
timestep 0.001
# 2nd equilibration run
compute ke all ke/atom
variable temp atom c_ke/(1.5*0.0000861) #Temperature in Kelvin for metal units, kb = 8.617e-5
fix 1 all nve
reset_timestep 0
compute layers all chunk/atom bin/1d x lower 1.824261252 units box
fix 2 all ave/chunk 1 1000 1000 layers v_temp file profilelong.mp
fix 3 all thermal/conductivity 5 x 20
variable tdiff equal f_2[11][3]-f_2[1][3]
thermo_style custom step temp epair etotal f_3 v_tdiff
thermo 1000
run 500000
# thermal conductivity calculation
# reset fix thermal/conductivity to zero energy accumulation
fix 3 all thermal/conductivity 5 x 20
fix ave all ave/time 1000 1 1000 v_tdiff ave running
thermo_style custom step temp epair etotal f_3 v_tdiff f_ave
run 15000000
***This gives the outputs:
f_3 = 1092964.2 eV
f_ave = average temperature difference = 512.12 K
total time of run = 15000 ps = 15ns
Dimension used: x = 36 A , y = 36 A , z = 84 A
Calculation:
dQ = 1092964.2 * 1.602e-19 / 15ns / 36 / 84 /(Angstroms)^2 / 2
dX = 36/2 Angstroms
dTemp = 512.12 K
Kappa = dQ x dX/dTemp ****
I should mention I did use the fix thermal conductivity command with No. of swaps equal to 10 as well, and the results that I obtained weren’t much different. I would be really grateful with some input on this method that I am following. Please let me know if I can provide more information.
Best regards,
Sayeed Nafis Sami.