Dear lammps-users,
I was experimenting around with the ‘in.langevin’ example file (/example/KAPPA folder) which explains how to calculate heat flux in a material block which has two sections maintained at different temperatures.
I have modified the code to calculate heat flux in a metal block, with two sections maintained at temperatures 300 K and 30 K.
I have attached the modified code below. The only modifications I have made are:
-
units have been changed to metals;
-
${rho} has been replaced with lattice constant a;
-
the variables t, {tlo} and {thi} have been changed to reflect the new temperatures;
-
pair_style and pair_coeff have been replaced with eam potential;
-
a minimization line has been included
-
the number of time steps in the run command to ensure equilibration
When I run the code, I notice that during the second equilibration, the values of {c_Thot} and {c_Tcold}, which I expected to equilibrate around 300K and 30 K, instead equilibrates around 190K and 125K. I have attached the relevant section (second equilibration) of the log file below:
Step c_Thot c_Tcold f_hot f_cold v_tdiff
10001 141.24863 150.10717 -0 -0 -8.8585426
11000 173.05511 120.18834 -11.692249 10.475999 52.866773
12000 193.60341 123.76728 -23.07725 20.586381 69.836134
13000 194.49871 121.18845 -33.856803 30.423574 73.310263
14000 193.43548 115.85178 -43.625538 39.474663 77.5837
15000 183.81897 121.83547 -54.938498 48.470818 61.983496
16000 187.36465 132.88362 -66.600173 58.805771 54.481023
17000 188.62838 122.57918 -77.151449 68.418573 66.049201
18000 196.42347 123.57975 -88.512999 78.067952 72.843719
19000 180.84652 123.11601 -97.683471 88.680379 57.730506
20000 187.91853 126.01149 -107.75662 98.596238 61.907041
21000 195.1383 135.0109 -117.57124 108.36986 60.127399
22000 206.50158 126.70738 -128.97982 118.37991 79.7942
23000 194.84583 128.95121 -141.14816 128.77797 65.894622
24000 190.39876 134.10713 -151.77277 138.30786 56.291625
25000 192.74748 137.06683 -162.96241 148.95368 55.680649
26000 206.78701 127.12741 -172.9532 159.59619 79.659593
27000 193.79335 121.94551 -183.33201 169.70739 71.847839
28000 194.8372 127.03072 -194.41835 179.87957 67.806484
29000 185.03134 124.83074 -205.17285 190.05637 60.200595
30000 187.42739 128.21589 -217.34754 200.45678 59.211508
30001 187.26346 127.67066 -217.33961 200.47427 59.592801
To fix this, what I did is to reduce the damping factor by a factor of 100 (from 1.0 to 0.01) - which by trial and error I found to be sufficient to get the temperatures to 300K and 30K as seen below.
Step c_Thot c_Tcold f_hot f_cold v_tdiff
10001 141.24863 150.10717 -0 -0 -8.8585426
11000 289.65325 33.44922 -58.660593 40.638714 256.20403
12000 294.67858 33.723322 -94.951599 69.17536 260.95526
13000 294.13983 32.105748 -126.39369 100.19958 262.03408
14000 284.19644 31.432519 -155.55538 127.36729 252.76392
15000 299.28698 35.014069 -182.14423 151.87195 264.27291
16000 301.87708 29.92551 -207.6585 177.17011 271.95157
17000 313.6155 33.249354 -232.53255 201.3258 280.36614
18000 294.57497 32.030982 -257.3655 224.18298 262.54399
19000 320.79769 31.339274 -283.24669 248.56029 289.45841
20000 303.50405 32.655081 -302.10371 273.17846 270.84897
21000 313.08154 32.027635 -328.58578 297.16818 281.0539
22000 302.36409 31.537685 -349.47573 319.77066 270.8264
23000 301.25575 30.48889 -374.34958 343.93773 270.76686
24000 287.28133 32.857736 -396.68116 368.3599 254.4236
25000 294.57522 30.480307 -424.40023 392.05486 264.09492
26000 295.73163 30.747492 -445.92789 415.63427 264.98413
27000 289.06742 33.475331 -468.67533 439.67745 255.59209
28000 308.06146 35.067424 -495.14317 462.55926 272.99403
29000 326.63032 33.976095 -518.30047 486.40306 292.65422
30000 296.69241 34.593652 -539.34018 510.14109 262.09876
30001 294.9412 33.997413 -539.19873 510.20073 260.94379
Q1) Is it fine to reduce the damping factor in fix langevin to any small value as needed in order to reach the temperatures that we require ?
Q2) I expected that the value of damping factor will not play a role in the final calculation of heat flux (as explained in the Readme file in the same folder). However, if I reduce the damping factor further from 0.01 to 0.005, the heat flux calculation reduces by a factor 1.7. How should I take into account the damping factor in the calculation of heat flux ?
The code I am using is:
sample LAMMPS input script for thermal conductivity of solid
thermostatting 2 regions via fix langevin
settings
variable x equal 10
variable y equal 10
variable z equal 20
variable a equal 3.52
variable t equal 150.0
variable rc equal 2.5
variable tlo equal 30.0
variable thi equal 300.0
setup problem
units metal
atom_style atomic
lattice fcc ${a}
region box block 0 $x 0 $y 0 $z
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create $t 87287
pair_style eam ## set interatomic potential style to be EAM
pair_coeff * * Ni_u3.eam ## read in interatomic potential file
neighbor 0.3 bin
neigh_modify delay 0 every 1
heat layers
region hot block INF INF INF INF 0 1
region cold block INF INF INF INF 10 11
compute Thot all temp/region hot
compute Tcold all temp/region cold
Minimization
minimize 1e-25 1e-25 5000 10000
1st equilibration run
fix 1 all nvt temp $t $t 0.5
thermo 1000
run 10000
velocity all scale $t
unfix 1
2nd equilibration run
fix 1 all nve
fix hot all langevin {thi} {thi} 1.0 59804 tally yes
fix cold all langevin {tlo} {tlo} 1.0 287859 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
variable tdiff equal c_Thot-c_Tcold
thermo_style custom step c_Thot c_Tcold f_hot f_cold v_tdiff
dump 2 all xyz 1000 dump2.xyz
thermo 1000
run 20000
thermal conductivity calculation
reset langevin thermostats to zero energy accumulation
compute ke all ke/atom
variable temp atom c_ke/1.5
fix hot all langevin {thi} {thi} 1.0 59804 tally yes
fix cold all langevin {tlo} {tlo} 1.0 287859 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
fix ave all ave/time 10 100 1000 v_tdiff ave running
thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave
compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.langevin
run 20000
Thanks!
Karthik