converting units from "real" to "metal"

Hello,

I have a question about converting units from "real" to "metal".
I've read similar kind of posts but I'm still wondering.

I converted units from "real" to "metal" in the sample input file in the
document of "compute heat/flux".

I use the same random number seed, and at the time step 0, the
calculated total energy is the same.
However, at timestep 2000, they have different temperature,
74.682256[K] from the "real" units input file and
66.078762[K] from the "metal" units input file.

and the calculated thermal conductivity have different result,
0.29[W/mK] from the "real" units input file and
0.19[W/mK] from the "metal" units input file.

Where are the differences from?
Anyone has any ideas?

My conversion process is below.
1. I converted values which have units of time[fs] to [ps] in the
definition of "dt" and "fs2s" and the Tdamp value in fix nvt.

2. I converted values which have units of energy[kCal/mol] to [eV] in
the definision of "kCal2J" and the 1st potential parameter in "pair_coeff".

3. I didn't change values which have units of temperature[K],
distance[Angstroms] and mass[g/mol] because "metal" and "real" has the
same units of them.

I will attach the both input file of "real" units and that of "metal"
units I used. (The run time is shortened from the original sample file)

#### converted input in "metal" units from here ####
# Sample LAMMPS input script for thermal conductivity of solid Ar

units metal
variable T equal 70
variable V equal vol
variable dt equal 0.004
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval

# convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.6021764E-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal \{eV2J\}\*{eV2J}/\{ps2s\}/{A2m}

# setup problem

dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.0103299 3.405
timestep ${dt}
thermo $d

# equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 0.01 drag 0.2
run 2000

# thermal conductivity calculation, switch to NVE if desired

#unfix NVT
#fix NVE all nve

reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix JJ all ave/correlate $s $p d &              c\_flux\[1\] c\_flux\[2\] c\_flux\[3\] type auto file J0Jt\.dat ave running variable scale equal {convert}/${kB}/$T/$T/$V*s\*{dt}
variable k11 equal trap(f_JJ[3])*\{scale\} variable k22 equal trap\(f\_JJ\[4\]\)\*{scale}
variable k33 equal trap(f_JJ[5])*${scale}
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run 10000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print "average conductivity: $k[W/mK] @ T K, {ndens} /A^3"
#### converted input in "metal" units till here ####

#### input in "real" units from here ####
# Sample LAMMPS input script for thermal conductivity of solid Ar

units real
variable T equal 70
variable V equal vol
variable dt equal 4.0
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval

# convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal \{kCal2J\}\*{kCal2J}/\{fs2s\}/{A2m}

# setup problem

dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.2381 3.405
timestep ${dt}
thermo $d

# equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 2000

# thermal conductivity calculation, switch to NVE if desired

#unfix NVT
#fix NVE all nve

reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix JJ all ave/correlate $s $p d &              c\_flux\[1\] c\_flux\[2\] c\_flux\[3\] type auto file J0Jt\.dat ave running variable scale equal {convert}/${kB}/$T/$T/$V*s\*{dt}
variable k11 equal trap(f_JJ[3])*\{scale\} variable k22 equal trap\(f\_JJ\[4\]\)\*{scale}
variable k33 equal trap(f_JJ[5])*${scale}
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run 10000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print "average conductivity: $k[W/mK] @ T K, {ndens} /A^3"
#### input in "real" units till here ####

Thank you in advance.

TADA Masahiro

I think this difference of physical properties between the results calculated from "metal units" and "real units" is not just a problem about thermal conductivity but about MD if my conversion is correct.

I am sorry I did not attach the input files before.
I attached input files and the outputs in this mail.
The included jpg image is about comparison of output of temperature, total energy, pair energy and pressure between metal units case and real units case.

In this results, the values are almost the same in the first 1000 steps. But temperature and total energy change the tendency with increase of time step.

Is this from summation of errors?
I considered Tdamp has time units, so I changed this parameter from fs to ps.

I show the differences of input files below.

$ diff in.Ar_metal.txt in.Ar_real.txt
##### differences from here #####
< units real

in.Ar_real.txt (1.97 KB)

properties_from_metal_units.txt (14.8 KB)

properties_from_metal_units_converted.txt (7.09 KB)

properties_from_real_units.txt (14.8 KB)

figs_unit_metal.jpg

in.Ar_metal.txt (1.96 KB)

If you do the same problem in 2 different units,
and if you have done all the conversions correctly,
then you should get the exact same answer, but
only for the first few 100 steps. After that
the 2 simulations will drift apart. That is the
nature of MD and integrating ODEs.

So what you are seeing looks normal.

Steve

Thank you very much for your reply.

I agree with you, it is from integrating ODEs.

I calculated in longer runtime with both input files on metal and real units and the results of thermal conductivity will converge the same value. This means they calculate the same MD simulations.

steps of run metal[W/mK] real[W/mK]
     10000 0.40 0.28
     20000 0.24 0.27
     30000 0.32 0.26
     50000 0.29 0.25

Thank you very much again.

TADA Masahiro

(2012/09/17 22:38), Steve Plimpton wrote:

I am sorry for that I mentioned only about the heat flux.
Other properties converged on almost the same values as well.

This is averages and standerd deviations of each property from 0 step to 250000steps.

                     average deviation
Temp[K] metal 69.87 3.33
            real 70.06 3.36
E_pair[J] metal -3.18E-18 3.44E-20
            real -3.18E-18 3.48E-20
TotEng[J] metal -2.81E-18 3.74E-20
            real -2.81E-18 4.02E-20
Press[Pa] metal 5.31E+07 1.89E+07
            real 5.57E+07 1.93E+07

These results mean they calculate the same MD simulations.

Thank you very much for your advice, Steve.

TADA Masahiro

(2012/09/18 4:29), TADA Masahiro wrote: