How to use fix thermal/conductivity and rNEMD method?

LAMMPS
1

Hello everyone,

I am trying to calculate thermal conductivity, using rNEMD method (Muller Plather algorithm) fix thermal/conductivity that is explained in LAMMPS manual.
To do so, after equilibrating the system, I ran NPT ensemble to have a stable system at 300 k. Then, I calculated kinetic energy, followed by fix NVE to get the production run and thermal conductivity.

My question is that does the procedure that I am following is correct? My concern mostly is for the fix NVE part, as I am not sure that for rNEMD method, should NVE ensemble be run or other ensembles, like NVT should be used to impose a temperature gradient to the system and then use fix thermal/conductivity along with that temperature gradient. Since, using NVE, I am not sure how the temperature gradient would be imposed to the system.

Also, I am using real units. So, is the calculation of kinetic energy correct? c_ke/0.00298?

Here is my input file that I am trying to run with NVE.

echo               both
units              real
dimension          3
boundary           p p p
atom_style         full

restart 50       restart

# Pair_Coefficient ..............................................................................................................................................


pair_style    	    lj/class2/coul/long 10 15
dielectric          1.0
bond_style          class2
angle_style         class2
dihedral_style      class2
improper_style      class2
#read_data          System.data
read_restart          System_2.data
#info coeffs out log

#...................................................................................................................................................................................

kspace_style        pppm 0.0001


neighbor           2  bin
neigh_modify       every 1 delay 0 check yes

#velocity        all create 300.0 497 dist gaussian


#minimize 1.0e-6 1.0e-6 5000 5000



timestep           1
thermo	           1000
#thermo_style       custom step cpu temp pe ke epair ecoul  etotal density lx ly lz

# reset timestep
reset_timestep 0

# saving trajectories
dump 1             all xyz 1000 System_2.xyz

#Exports ......................................................................................................................................

#nvt run
#fix 2 all          nvt temp 300.0 300.0 50.0
#run	               400000
#undump
#unfix 2

#..............................................................................................................................................
#raise temperature
#fix 3             all nvt temp 300.0 600 50
#write_restart 	   restart-1.equil
#run               1000000
#undump            1
#unfix 3
#..............................................................................................................................................

#run at high temperature for 3000ps
#fix 4              all npt temp 600 600 50 iso 1 1 1000
#write_restart 	   restart-1.equil
#fix nvt all nvt temp 600.0 600 50
#run                3000000 upto
#undump             1
#unfix 4
#..............................................................................................................................................

# anneal
#fix 5 all nvt temp 600.0 300 50
#run 300000  #ADJUST the steps here for different cooling rates
#unfix 5

#..............................................................................................................................................

#fix 6 all npt temp 300 300 50 iso 1 1 1000
#fix nvt2 all nvt temp 570.0 570.0 50
#run 3000000 upto #ADJUST STEPS of production run
#unfix 6
#write_data Nafion_3landa.data

#..............................................................................................................................................

compute	ke all     ke/atom
variable temp atom c_ke/0.00298  


fix 7 all         nve

#calculate temperature difference
compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 8 all          ave/chunk 100 10 1000 layers v_temp file profile.mp

# Calculate thermal conductivity;   velocity exchange is 10 here
fix		9 all        thermal/conductivity 10 z 20
variable           tdiff equal f_8[11][3]-f_8[1][3]

# thermal conductivity calculation
# reset fix thermal/conductivity to zero energy accumulation

fix		10 all        thermal/conductivity 10 z 20
fix   ave all      ave/time 1 1 1000 v_tdiff ave running
thermo_style       custom step cpu temp pe epair etotal density f_10 v_tdiff f_ave
dump 1             all xyz 1000 System3.xyz
run	               1000000
undump             1
write_restart 	   restart-3.equil

I really appreciate your guidance on this matter as I could not find much information in this regard.

Cheers,

Mahsheed

2

In LAMMPS those are not “ensembles” but “integrator fix commands”. In fact, neither will produce the indicated ensemble, because those require a system in equilibrium.

What is correct should be looked up in the publication(s) describing the method. You can use the examples that LAMMPS provides as guidelines, but neither the manual nor the examples are a replacement for properly studying the publications of the method (and possibly related text books if available).

3

Dear Mahseed,
Yes, you need to use the NVE ensemble (although, as Dr. Kohlmeyer mentioned, they have a different meaning in LAMMPS), as the idea of the Muller-Plather method is to keep the energy of the system constant (loss of energy would result in error). I highly recommend you read their paper on the matter (cited in lammps docs), although it looks a bit complicated at the beginning, it makes much sense after you understand it.
In addition, you need to read recent articles on the topic as there are several considerations while using this method to achieve reasonable results (You can read my article on the subject, but don’t limit yourself to it and read other works as well.)

1 Like
4

Dear @anazarahari , thank you for your reply. I checked the paper and also I am reading your interesting paper. As you said, the idea behind this technique now is making more sense to me.

However, for the calculation part, I am not sure of the conversion factor. I am using the same technique as described in here,but I am not sure if this is a correct way to convert the obtained TC from real units to W/(m.k). Could you also comment on this part if possible for you?

Thanks,

Mahsheed

5

Hi Again,
I am not familiar with the unit conversions of “real” units. However, I believe it doesn’t make much difference as long as you know the unit of the results that LAMMPS produces. You can search on the internet and find the conversion values.

Saying this, As far as I remember, you need to convert time (to seconds), energy (to joles), and length (to meters). Remember that if you are using the LAMMPS example file for thermal conductivity, the file was prepared for an LJ system. So you need to make the necessary adjustments for that as well.