Question about viscosity calculation using R-NEMD method

Dear Lammps developers,

I am currently using R-NEMD methods for viscosity calculation.

According to some exited papers, viscosity of Cu50Zr45Al5 at T=900K was above 100 CP. However, my calculation result of Cu45Zr45Al10 was only 30 CP.

I’ve made a lot of tests, but seems it’s hard to achieve a high value of viscosity at 900K and a linear profile at the same time.

I have tested a series of parameters in “fix viscosity” command to adjust the frequency of momentum exchanges (Nevery and Nswap).

When I set Nevery =1 and Nswap =2, the value of viscosity was 30 CP, and the velocity profile was almost linear.

However, when I set Nevery =1 and Nswap =1, the value of viscosity became 78 CP, but the velocity profile was no longer linear.

It was the same when I set Nevery =100 and Nswap between 50 and 100, and the value were 158 CP and 58 CP, respectively. But the former was highly non-linear.

Could you give me some advice about the above situations?

Besides, here are some efforts I’ve made to ensure the correctness:

(1) I’ve tested the influence of cooling rate, relaxation time, potential, simulation box size and other parameters, but none of them seems to matters.

(2) The values of viscosity at temperatures from 900 to 1600 K can be well described by the VFT equation.

(3) I have also conducted the simulation using the Green-Kubo methods, and the value at T=900K was between 20 to 30 CP and the values at other temperatures were also close to R-NEMD’s.

Thank you.

Sincerely,

Xiaochang Tang

My calculation files:

  1. in.mp.3d

read_restart restart.900000

variable t equal 900

------------------------ FORCE FIELDS ------------------------------

pair_style eam/alloy

pair_coeff * * ZrCuAl(2011).eam.alloy Zr Cu Al

mass 1 91.224

mass 2 63.546

mass 3 26.982

------------------------ EQUILIBRATION ------------------------------

reset_timestep 0

timestep 0.001

fix 1 all nvt temp $t $t 0.1

thermo 1000

thermo_style custom step temp ke pe etotal press

run 10000

------------------------ Muller-Plathe ------------------------------

fix 4 all viscosity 1 x y 20 swap 2

compute apple all chunk/atom bin/1d y center 0.05 units reduced

fix 5 all ave/chunk 10 100 1000 apple vx density/number temp file profile.mp.3d.${t}K

variable dVx equal f_5[11][3]-f_5[1][3]

thermo 1000

thermo_style custom step temp ke pe etotal press f_4 v_dVx

run 10000

------------------------ Date gathering ------------------------------

fix 4 all viscosity 1 x y 20 swap 2

variable ps2s equal 1.0e-12

variable A2m equal 1.0e-10

variable g2kg equal 1.0e-3

variable Avogadro equal 6.022e23

variable convert equal ({g2kg}/{Avogadro})/({ps2s}*{A2m})

variable visc equal -((f_4/(2*(step0.001-20)lxlz+1.0e-10))/(v_dVx/(ly/2)))${convert}

fix vave all ave/time 1000 1 1000 v_visc ave running start 21000

thermo_style custom step temp ke pe etotal press f_4 v_dVx v_visc f_vave

variable vave equal f_vave

run 1000000

print “{t} {vave}” append viscosity.txt

  1. in.cooling

------------------------ INITIALIZATION ----------------------------

units metal

dimension 3

boundary p p p

atom_style atomic

----------------------- ATOM DEFINITION ----------------------------

lattice fcc 3.52

region a block 0 30 0 30 0 30

create_box 3 a

create_atoms 1 region a

------------------------ FORCE FIELDS------------------------------

pair_style eam/alloy

pair_coeff * * …/ZrCuAl(2011).eam.alloy Zr Cu Al

mass 1 91.224

mass 2 63.546

mass 3 26.982

set type 1 type/fraction 2 0.45 94756789

set type 1 type/fraction 3 0.181818182 87899075

------------------------ SETTING -----------------------------------------

timestep 0.001

compute RDF all rdf 500 1 1 1 2 1 3 2 2 2 3 3 3

thermo 10000

thermo_style custom step temp vol pe ke press enthalpy epair

variable temp equal temp

variable pe equal pe

variable vol equal vol

variable enthalpy equal enthalpy

------------------------ GLASS FORMING ------------------------------

velocity all create 1800 5812384 rot yes dist gaussian

fix 1 all npt temp 1800 1800 0.1 iso 0 0 1

run 500000

unfix 1

write_restart restart.1800K

reset_timestep 0

fix 1 all npt temp 1800 300 0.5 iso 0 0 1

fix rdf all ave/time 100000 1 100000 c_RDF file 1800-300k.rdf mode vector

fix Tg all print 10000 “{temp} {pe} {vol} {enthalpy}” file Tg.txt screen no

restart 100000 restart.*

dump 1 all custom 100000 dump.* id type x y z

dump_modify 1 element Zr Cu Al

dump_modify 1 sort id

run 1500000

print “All done!”