LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] tally fix langevin

# Re: [lammps-users] tally fix langevin

 From: Steve Plimpton Date: Thu, 12 Oct 2017 07:14:00 -0600

The doc page for fix langevin explains that
the damping factor has units of time.
So if you leave it at 1.0 (the LJ value),
that means your T will relax on the order
of a picosec (real units), which is a lot
of timesteps.  So ye, adjusting the damping
factor is normally required when changing
units.

The damping factor is not part of the heat
flux calculation, but that assumes you
have run long enough to be in equilibrium.

Steve

On Wed, Oct 4, 2017 at 5:30 AM, Karthik Sasihithlu wrote:
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

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most