We are currently trying to calculate the lattice thermal conductivity of several metals (Ni, Cu, Pd, Ag, Pt and Au) using the non-equilibrium Müller-Plathe method in LAMMPS. We are considering big supercells (20x20x20 unit cells, 70k+ atoms) for each elemental metal separately. The interactions between atoms are mediated by the well-known MEAM potential. The heat flux is tallied using the fix thermal/conductivity command, and temperature profiles are recorded using the compute chunk/atom and fix ave/chunk commands. Please find attached the typical input run for reference. We have successfully run this for other systems (molecular liquids governed by OPLS-AA), but for MEAM metals we found that the temperature profiles are basically erratic noise. While the lattice thermal conductivity of metals should low, we do not expect it to be completely negligible.
Is anyone aware of existing problems with the application of the Müller-Plathe method together with the MEAM potential? Is there particular requirements for such a case? We have played with the Nevery and Nswap parameters in the fix thermal/conductivity command with no success.
Thank you in advance for your attention.
*****
include "system.in.init"
read_data "system.data"
include "system.in.settings"
include "system.in.charges"
neigh_modify every 1 delay 0 check yes
variable t equal 0.001
variable T equal 298.0
variable P equal 1.013
variable s equal 5
variable c equal 10000
variable d equal $s*$c
variable kB equal 8.617333262E−5
minimize 1.0e-6 1.0e-8 1000 100000
reset_timestep 0
timestep $t
thermo_style custom step etotal temp press lx ly lz density
thermo 1000
fix 1 all nvt temp $T $T 10.0
run 500000
unfix 1
fix 2 all npt temp $T $T 10.0 iso $P $P 100.0
run 2000000
unfix 2
reset_timestep 0
thermo_style custom step etotal epair ke temp press
thermo 1000
fix 3 all nvt temp $T $T 10.0
run 10000000
unfix 3
reset_timestep 0
fix 4 all nve
fix 5 all thermal/conductivity 100 z 50
compute ke all ke/atom
variable temp atom c_ke/(1.5*${kB})
compute layers all chunk/atom bin/1d z lower 0.02 units reduced
fix MP all ave/chunk $s $c $d layers v_temp file temp.profile ave one
thermo_style custom step etotal epair ke temp press f_5
thermo 1000
run 10000000
unfix MP
unfix 5
unfix 4
I think for OPLS-AA, you are using 'real' unit and for MEAM it is 'metal'. Those are not equivalent! So, you have to change the tdamp and pdamp values in npt/nvt fixes.
Dear Tanvir Ahmed,
Thanks for your answer. We are using 'units metal' in our simulations. You can infere that from the variables we have defined above. It is true that the Tdamp and Pdamp values should be accordingly modified. By using Tdamp = 10 in 'units metal' (i.e. 10 ps) it takes much longer for the system to equilibrate temperature than the same Tdamp value in 'units real' (i.e. 10 fs) (which we have verified). However, the Muller-Plathe thermal conductivity calculations are done in NVE, so once equilibrated that should not affect.
Happy to hear if you have any other alternatives.
- Badges
- Science topic
- Similar topics
- Biological Science
- Molecular Cell Biology
- Culture Cells
- Cell Line