Regarding Diffusion Coefficients With Lammps

LAMMPS LAMMPS Mailing List Mirror
1

Hello All,

Maybe you all can help me solve my problem. I have been trying to get diffusion coefficients from my Cementite system through several methods. Using the Test file in the LAMMPS Example Directory ATC_Twotemperature I extracted the MSDs at multiple temperatures and various timescales. (Using real units) For long enough timescales there were no problems. However when I try to do the same on my cementite system I get RSMD that even with a log log graph doesn’t make much sense at any temperature. I am still testing timescales. (Metal units). I am posting a picture of what the data looks like over time. Any ideas as to why it might look like this?

importantsystem.png

2

four comments:

  1. MSD != RMSD
  2. the example system you are referring to is simulating liquid argon, which is very different from your material.
  3. when you visualize your system, is it liquid? i.e. do atoms actually move around?
  4. there is not enough information here to provide any more specific recommendations

Axel.

3

Hi Axel,

Yes I know that MSD =! RMSD. Still the slope from MSD does not give me a value that makes sense.
Yes I know that system is argon. But I want to calculate diffusion coefficient after the effect of using ATC so i wanted to see if the process works
Yes I visualized the system. No Its not liquid yes but the atoms are moving. It should be liquid above 1700 K

I have attached my file here.

#Making Melt Fe3c Modeled after in.elastic and Melt File made by Carlos Ruestes

units metal
atom_style full
boundary p p p
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no

read_data data.Fe3C.xtal.3x3x3
change_box all triclinic

pair_style meam/c
pair_coeff * * Fe3C_library_Liyanage_2014.meam Fe C Fe3C_Liyanage_2014.meam Fe C

Define minimization parameters

variable etol equal 1.0e-3
variable ftol equal 1.0e-10
variable maxiter equal 100
variable maxeval equal 1000
variable dmax equal 1.0e-2

Define MD parameters

variable nevery equal 10 # sampling interval
variable nrepeat equal 10 # number of samples
variable nfreq equal {nevery}*{nrepeat} # length of one average
variable nthermo equal {nfreq} # interval for thermo output variable nequil equal 10*{nthermo} # length of equilibration run
variable nrun equal 3*${nthermo} # length of equilibrated run
variable inittemp equal 1200.0 # temperature of initial sample
variable hightemp equal 1200.0 # temperature of above melting sample
variable timestep equal 0.001 # timestep
variable mass1 equal 28.06 # mass
variable adiabatic equal 0 # adiabatic (1) or isothermal (2)
variable tdamp equal 0.01 # time constant for thermostat
variable seed equal 123457 # seed for thermostat
variable tdamp equal “v_timestep10" # May require some systematic tests for proper adjustment
variable pdamp equal "v_timestep100” # May require some systematic tests for proper adjustment

Setup minimization style

min_style cg
min_modify dmax ${dmax} line quadratic

thermo 100
thermo_style custom step temp vol press ke pe etotal

Compute initial state

fix 3 all box/relax aniso 0.0
minimize {etol} {ftol} {maxiter} {maxeval}
velocity all create {inittemp} {seed} rot yes dist gaussian

dump id all atom 500 dump.melt5

reset_timestep 0

fix equilibration all npt temp {inittemp} {hightemp} {tdamp} iso 0 0 {pdamp}
run 2000000

#fix equilibration all npt temp {hightemp} {inittemp} {tdamp} iso 0 0 {pdamp}
#run 20000

#fix equilibration all npt temp {inittemp} {inittemp} {tdamp} iso 0 0 {pdamp}
#run 20000

write_data data.meltFe3C.* # data file for restarts

print “all DONE”
write_restart meltFe3c.r

4

if it is not liquid how should there be significant diffusion?

5

I dont think there should be significant diffusion. But something measurable. I will retry at even higher temperatures and get back to you.

6

Hey Axel,

Ok I ran this at a higher temperature and I’m still getting this graph that just doesn’t make sense to me.

7

Since you have yourself confirmed that the software itself is working correctly (and I have myself as well as have many others), the most likely conclusions are that either a) your system is still not liquid and thus not diffusing as expected by the Einstein relation that connects the MSD to the self-diffusion coefficient (technically speaking the derivation was made for brownian motion of gases) or b) there is something incorrect with your force field parameters or simulation settings (e.g. incorrect units or unit conversions).

Axel.

8

You can always do a sanity test yourself by dumping 2
snapshots of atom coords at time 0 and some later time
(with image flags in the snapshots), and computing
the MSD yourself as a post-processing operation. Then
compare that single value to what compute msd is giving
you at the same point in time. If you viz the 2 snapshots
you can likely also understand how much atoms are
moving (or not).

Steve