LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] tally fix langevin
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] tally fix langevin


From: Karthik Sasihithlu <skarthik.rao@...24...>
Date: Wed, 4 Oct 2017 12:30:57 +0100

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